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^  In  a  recent  paper  (see  [11]),  L.  Greengard  and  V.  Rokhlin  introduce  a  numerical  technique  for 
the  rapid  solution  of  integral  equations  resulting  from  linear  two-point  boundary  value  problems 
for  second  order  ordinary  differentia)'  ^nations.  In  this  paper,  we  extend  the  method  to  systems 
of  ordinary  differential  equations.  Aftfer  reducing  the  system  of  differential  equations  to  a  system 
of  second  kind  integral  equations/  discretize  the  latter  via  a  high  order  Nystrom  scheme.  A 
somewhat  involved  analytical  appai/atus  is  then  constructed  which  allows  for  the  solution  of  the 
discrete  system  using  0{N  •  ■  n^)  operations,  with  N  the  number  of  nodes  on  the  interval, 

p  the  desired  order  of  convergence,  and  n  the  number  of  equations  in  the  system.  Thus,  tlie 
advantages  of  the  integral  equation  formulation  (small  condition  number,  insensitivity  to  boundary 
layers,  insensitivity  to  end-point  singularities,  etc.)  are  retained,  while  achieving  a  computational 
efficiency  previously  available  only  to  finite  difference  or  finite  element  methods. 

We  in  addition  present  a  Newton  method  for  solving  boundary  value  problems  for  nonlinear  first 
order  systems  in  which  each  Newton  iterate  is  the  solution  of  a  second  kind  integral  equatio.r,  the 
analytical  and  numerical  advantages  of  integral  equations  are  thus  obtained  for  nonlinear  boundary 
value  problems. 
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I.  Introduction 


Second  kind  integral  equations  have  been  a  popular  analytical  tool  in  the  study  of  ordi¬ 
nary  differential  equations  for  nearly  a  century.  When  boundary  value  problems  are  being 
considered,  the  integral  equations  which  arise  are  of  the  Fredholm  type.  From  an  abstract 
viewpoint,  the  advantage  of  this  formulation  is  that  many  properties  of  the  solution  are  readily 
apparent;  from  a  computational  viewpoint,  the  linear  systems  which  arise  from  discretization 
are  generally  well-conditioned.  An  ill-behaved  differential  equation  can  often  be  reduced  to  a 
perfectly  tractable  integral  equation  by  means  of  an  appropriate  choice  of  the  “background” 

Green’s  function.  Standard  finite  difference  and  finite  element  methods,  on  the  other  hand, 
which  discretize  the  original  differential  equation,  encounter  serious  numerical  difficulties  when 
the  solution  possesses  derivatives  of  large  magnitude  (boundary  layers).  A  second  advantage  is 
that  there  exist  extremely  stable,  high  order  numerical  methods  for  the  solution  of  second  kind 
Fredholm  equations,  while  the  order  of  convergence  of  most  practical  schemes  for  the  solution 
of  ordinary  differential  equations  tends  to  be  limited,  even  if  Richardson  extrapolation  and 
deferred  correction  approaches  are  considered. 

Despite  aU  these  advantages,  integral  equations  are  virtually  never  used  as  a  numerical  tool 
for  the  solution  of  two-point  boundary  value  problems,  since  their  discretization  leads  to  dense 
systems  of  linear  algebraic  equations,  and  the  solution  of  a  dense  linear  system  of  dimension  A' 
requires  order  0{N^)  arithmetic  operations.  Finite  difference  and  finite  element  schemes  lead 
to  banded  systems  of  linear  algebraic  equations,  and  the  solution  of  the  latter  requires  order 
0{N)  arithmetic  operations,  with  N  the  dimension  of  the  problem.  This  makes  the  use  of 
integral  equations  extremely  unattractive  as  a  numerical  tool,  despite  their  superior  analytical 
properties.  A  similar  difficulty  is  encountered  when  spectral  methods  are  applied  to  boundary 
value  problems.  They  yield  high  order  accuracy,  but  result  in  dense  systems  of  linear  algebraic 
equations. 

Recently,  [11]  presented  a  fast  numerical  algorithm  for  solving  two-point  boundary  value 
problems  for  second  order  differential  equations.  By  solving  the  problems  as  second  kind 
integral  equations,  one  obtains  the  superior  properties  of  integral  equations  over  differential 
equations.  By  using  the  technique  of  [11],  integral  equations  arising  from  boundary  value 
problems  are  solved  in  order  0(N)  arithmetic  operations. 

In  this  paper,  we  extend  the  results  of  [ll]  by  showing  that  integral  equations  arising  from 
two-point  boundary  value  problems  for  systems  of  ordinary  differential  equations  can  also  be 
solved  for  a  cost  proportional  to  the  number  of  nodes  N.  We  in  addition  present  a  Newton 
method  for  solving  boundary  value  problems  for  nonlinear  first  order  systems  in  which  each 
New'ton  iterate  is  the  solution  of  a  second  kind  integral  equation. 

The  plan  of  this  paper  is  as  follows;  in  Section  2  we  summarize  both  the  theory  of  Green's 
functions  for  first  order  linear  systems  and  the  theory  of  Newton  methods  for  first  order  non¬ 
linear  systems,  in  Section  3  we  develop  the  analytical  apparatus  to  be  used,  and  in  Section  4 
we  describe  the  numerical  schemes  themselves.  The  performance  of  the  methods  is  illustrated 
in  Section  5  with  numerical  examples.  Our  conclu.sions  are  discussed  in  Section  6. 

The  present  paper  is  similar  to  [11]  in  that  while  it  is  based  on  a  sequence  of  fairly  simple 
observations,  the  details  of  the  algorithm  are  somewhat  involved.  We  attempt  in  this  paper  to 
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present  both  cursory,  qualitative  descriptions  as  well  as  detailed,  rigorous  proofs. 


II.  Mathematical  Preliminaries 

In  this  section,  we  summarize  the  relevant  properties  of  both  the  boundary  value  problems 
to  be  addressed  and  the  second  kind  integral  equations  to  be  used  for  their  solution.  Most  of  the 
results  are  classical  and  can  be  found,  for  example,  in  [3]  and  [4].  The  rest  are  straightforward 
generalizations  to  systems  of  ordinary  differential  equations  of  well-known  facts  concerning 
second  order  boundary  value  problems  (see,  for  example,  [5]). 

2.1.  Notation  and  Definitions 

Definition  2.1  A  linear  first  order  system  of  ordinary  differential  equations  is  an  expression 
of  the  form 

$'(i)-|-p(i)-$(.t)  = /,  (1) 

with  $  :  [a,c]  R"  in  C^[a,c],  p  :  [a,c]  —  L(R”^”)  and  f  :  [a,c]  ^  R"  continuous,  and 
j^(RnxTi)  ihg  linear  space  of  all  linear  operators  R"  — >  R". 

Definition  2.2  If  f{x)  =  0,  (1)  assumes  the  form 

$'(i)  +  p(i)  •  $(x)  =  0,  (2) 

and  is  referred  to  as  a  linear  homogeneous  first  order  system  of  ordinary  differential  equations. 

Definition  2.3  A  differentiable  function  $  ;  [a,c]  — ♦  R”  is  a  solution  to  a  linear  first  order 
boundary  value  problem  if  it  satisfies  an  equation  of  the  form  (1),  subject  to  boundary  conditions 
of  the  form 

/I  •  $(a) -h  C  •  4>(c)  =  7.  (3) 

with  A,C  £  L(R"’^”),  and  7  G  R". 

Definition  2.4  7/7  =  0,  (3)  becomes 


A  -^(a)  +  C  -^{0)  =  0.  (4) 

and  is  referred  to  as  a  set  of  homogeneous  boundary  conditions. 

Definitions  2. 5-2. 6  are  the  nonlinear  analogues  to  Definitions  2.1  and  2.3. 

Definition  2.5  A  nonlinear  first  order  system  is  defined  as  an  expression 

4>'(x)  =  F($(x),  .t),  (3) 

with  $  ;  [a,c]  — ►  R”  in  C^la,c],  F  :  R""^*  — ►  R”  continuous. 
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Definition  2.6  A  differentiable  function  $  :  [a,c]  — ♦  R"  is  a  solution  to  a  nonlinear  first 
order  boundary  value  problem  if  it  satisfies  an  equation  of  the  form  (5),  subject  to  boundary 
conditions  of  the  form 

A  •  $(a)  +  C  •  $(c)  =  7,  (C) 

with  A,C  6  L(R”''"),  7  e  R". 

Definition  2.7  A  continuous  function  G{x,t)  :  [a,c]x{a,c]  — *  L(R"’^”)  is  the  Green  s  function 
for  a  boundary  value  problem  (1),  (4)  if 
1-  is  continuous  except  at  x  =  t, 

2.  G{x  +  0,i)  -  G{x  -  0,i)  =  In  for  all  x  £  [a,c], 

3.  ^G{x,t)  +  p(ar)  •  G{x,t)  =  0  for  all  x,t  £  [a,c],2:  #  t, 

4.  A  ■  G(a,i)  +  C  ■  G(c,t)  =  0  for  all  t  £  [a,c]. 

Remark  2.1  Green’s  functions  are  the  principal  analytical  tools  which  enable  boundary  value 
problems  to  be  solved  as  second  kind  integral  equations.  However,  Green’s  functions  are  known 
or  computable  for  very  few  problems  (1),  (4).  Fortunately,  we  can  use  one  of  the  known  Green  s 
functions  when  constructing  the  second  kind  integral  equation  for  a  particular  boundary  value 
problem.  When  a  Green’s  function  unrelated  to  a  problem  (1),  (4)  is  used  to  convert  that 
problem  to  an  integral  equation,  we  will  refer  to  this  Green’s  function  as  a  background  Green  s 
function. 

Definition  2.8  A  function  T  :  [a,c]  -♦  L(R”*'’’)  is  called  a  fundamental  matrix  for  (2)  if  it  is 
nonsingular  and 

T'(i)  +  p(a:)-T(x)  =  0  (T) 

for  all  X  £  [a,  c]. 

We  define  boundary  condition  matrices  D  and  to  be  used  in  theorems  in  the  remainder 
of  Section  11. 

Definition  2.9  Given  a  fundamental  solution  matrix  T  of  the  system  (2),  and  a  pair  of  ma¬ 
trices  A,C  given  by  (3),  the  boundary  condition  matrices  D,Dn  G  L(R”^’’)  are  defined  by  the 
formulae 

n  =  A  ■r(a)  +  C -Tic),  (S) 

L>^  =  A  +  C.  (9> 

We  define  a  residual  mapping  A’  and  Newton  iterates  to  be  used  in  a  Newton  method 
for  nonlinear  boundary  value  problems. 

Definition  2.10  Given  functions  Go,Pq  :  [a,c]  X  {a,c]  — ♦  L{R’”'’‘),  we  define  the  residual 
mapping  K  :  — >  R’‘  by  the  formula 

K{a{x),x)  =  a(x)  -  pn(x)  ■  f  GoixJ)  •  er(l)dt  -  F  (  f  Go(3:,t)  ■  cr^t)  d1 .  x^  .  (10) 
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Definition  2.11  For  any  continuous  ctq  :  [a,  c]  -+  R”,  we  refer  to  the  continuous  functions 
6ic  :  [a,  c]  — »  R"  as  Newton  iterates  if  for  each  k  =  1,2,..., 

(11) 

with  each  continuous  :  [a,  c]  —*  R”  recursively  defined  via  the  formula 

•  (g;c+i(j)  -  gfc(y))  =  -A'(<Tfc(i),i),  1:  =  0,1,---.  (12) 

dok 

Finally,  we  define  a  transposition  operator  to  be  used  in  the  paper. 

Definition  2.12  Given  an  interval  [61,62]  C  R  and  a  mapping  x  :  ^^[61,62]  — ►  L(R"’‘").  thf 
transpose  :  (^^[61,62])"  —  R'*  of  x  is  defined  by  the  formula 

X^('7)  =  ^  xit)-(Tii)dt,  (13) 

with  a  E  (A^)". 


2.2.  Green’s  Functions  for  First  Order  Systems 

Theorems  2. 1-2.8  provide  the  tools  for  the  conversion  of  first  order  systems  of  diffciential 
equations  into  second  kind  integral  equations.  Theorems  2.1,  2.2,  2.5  and  2.6  are  well  known 
and  can  be  found,  for  example,  in  [3]  and  [4].  The  authors  failed  to  locale  the  remaining 
theorems  in  the  literature. 

Theorems  2. 1-2.2  provide  conditions  for  the  existence  and  uniqueness  of  solutions  to  (1). 

(4). 

Theorem  2.1  For  any  continuous  function  p  :  [a,c]  — *  Z-(R”’'"),  the  homogeneous  first  order 
system  (2)  has  exactly  n  linearly  independent  solutions. 

Theorem  2.2  If  the  matrix  D  defined  by  (8)  is  nonsingular,  then  there  is  a  unique  solution 
to  the  equation  (1)  satisfying  homogeneous  boundary  conditions  (4).  Furthermore,  the  solutiem 
to  the  homogeneous  equation  (2)  satisfying  homogeneous  boundary  c^nd  ‘  ms  (4)  is  4>(.i  )  =  0. 

The  purpose  of  the  following  two  theorems  is  to  permit  the  conversion  of  problems  with 
inhomogeneous  boundary  conditions  to  those  with  homogeneous  ones.  Theorem  2.3  coll(el  Il^ 
linear  problems  of  the  form  (1),  (3);  Theorem  2.4  concerns  nonlinear  problems  of  the  form  (.5). 
(6). 

Theorem  2.3  If  the  boundary  condition  matrices  D,Dn  defined  by  (S),  (9)  are  both  nonsin¬ 
gular,  then  the  solution  to  the  problem  (1),  (S)  is  given  by  the  formula 

$(1)  =  $(i)  -f  I/,  (14) 
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with  V  €  R”  given  by  the  formula 


//  =  (A  +  C)-' .7,  (1-^) 

and  $  :  [a,c]  -+  R"  in  C^[a,c]  the  solution  to  the  first  order  system 

$'(1)  +  p(x)  •  $fj)  =  f{x)  -  p{x)  ■  u,  (IG) 

satisfying  homogeneous  boundary  conditions  (4). 

Proof.  Since  the  matrix  D  is  nonsingular,  it  immediately  follows  from  Theorem  2.2  that  there 
exists  a  unique  satisfying  the  equation  (16).  Substituting  (14)  into  boundary  conditions  (3), 
we  obtain 

v4  ■  (l'(a)  +  I/)  +  C  •  ($(c)  +  t-)  =  7.  (17) 

Now,  (15)  is  easily  obtained  from  the  combination  of  (17)  and  (4),  while  (16)  is  a  result  of 
substituting  (14)  into  (1).  □ 

Theorem  2.4  If  there  exists  a  unique  solution  ;  [a,c]  — ►  R”  to  the  problem  (5),  (6).  and  if 
the  matrix  D;^’  defined  by  (9)  is  nonsingular,  then  $  is  given  by  the  formula 

4>(i)  =  $(1)  +  t/,  ( IS) 

with  u  €  R"  given  by 

i/  =  (.4  +  C)-*.7,  (19) 

arid  ;  [fl,c]  — ’  R"  6  C'[a,c]  the  solution  to  the  nonlinear  boundary  value  problem 

4>'(x)  =  F(4>  +  i/,a:)  (20) 

with  homogeneous  boundary  conditions  (4). 

Proof.  Substituting  (18)  into  boundary  conditions  (6)  we  obtain 

>l-($(a)  +  i/)  +  C-(^(c)  +  i/)  =  7.  (21) 

Now.  (19)  is  easily  obtained  from  the  combination  of  (21)  and  (4),  while  (20)  is  a  result  of 
substituting  ( IS)  into  (5).  □ 

Theorem  2.5  provides  an  explicit  construction  for  the  Green’s  function  for  a  boundary  vahie 
problem  with  a  known  fundamental  matrix  T.  Given  a  Green’s  function  for  a  lioniogeneous 
problem  (2).  (4),  Theorem  2.6  provides  an  explicit  solution  for  the  inhomogeneous  problem  ( 1 ). 

(-I)- 

Theorem  2.5  If  the  matrix  D  defined  by  (8)  is  nonsingular,  then  there  exists  a  unique  Green's 
function  G  ;  [a,c]  X  [a,c]  — ♦  L(R"^”)  for  (2),  (4)-  G  is  given  by  the  formula 


-H 


T(i).(T-i(0  + J(0) 
r{x).j{t)  it>x), 
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with  J  :  [a,c]  — »  L(R"’‘")  given  by  the  formula 

J(0  =  -£>■' •C-T(c)-T-'(0,  (‘-^3) 

andr  :  [a,c]  ^  )  the  fundamental  matrix  for  (2)  (see  Definition  2.8). 

Theorem  2.6  Given  a  Green’s  function  for  the  problem  (2),  (4),  the  solution  $  for  the  problem 
(1),  (4)  can  be  obtained  via  the  formula 

$(i)=  r  G{x,t).fit)dt.  (24) 

Ja 

Tlie  following  two  theorems  are  the  principal  analytical  tools  used  in  this  paper.  Theorem 
2.7  is  used  to  reduce  a  linear  boundary  value  problem  (1),  (4)  to  a  second  kind  integral  equation, 
even  when  the  Green’s  function  for  the  problem  is  not  available;  Theorem  2.S  is  used  in  the 
same  fashion  to  reduce  nonlinear  boundary  value  problems  (5),  (4)  to  nonlinear  second  kind 
integral  equations. 

Theorem  2.7  Suppose  po  :  [a,c]  — ^  L(R"’‘")  is  continuous,  To  :  [o.c]  —  L(R’'‘’‘’‘)  tlu. 
fundamental  matrix  for  the  equation 

^'(x)  +  poix)-^(x)  =  0,  (25) 

and  Go  :  [a,  c]  X  [a,  c}  -*  L(R"’'”)  is  the  Green’s  function  for  the  boundary  value  problem  (25). 
(4).  Suppose  further  that  the  matrix  D  defined  by  (8)  and  the  matrix  Do  €  L(R’‘’‘")  defined 
by  the  formula 

Do  =  ^-To(a)  +  C-To(c),  (2G) 

are  both  nonsingular.  Then  the  solution  $  to  the  problem  (1),  (4)  can  be  obtained  via  the. 
formula 

$(i)=  r  Go{x,t)-a{i)dt,  (27) 

J  a 

with  a  :  [a.c]  — '  R"  the  solution  to  the  second  kind  integral  equation 

<rix)i-[]i{x)-pQ{x)}-  I  Go(x,t)  ■  a{t)dt  =  fix).  (2S) 

Ja 

Proof.  By  Theorem  2.2,  if  matrices  D,Do  are  nonsingular  then  the  problems  (1),  (4) 
and  (25),  (4)  have  unique  solutions,  and  therefore  the  background  Green’s  function  Go  i-''  also 
unique,  and  is  defined  by  Theorem  2.4.  Now,  (28)  is  obtained  by  substituting  (27)  into  (1).  □ 

Remark  2.2  Ifpo(^)  =  p(2:),then  the  solution  toequation  (28)  is  trivially  o  —  f .  Our  working 
assumption  is  that  for  some  background  problem  (25),  (4),  the  Green’s  function  is  known  or 
computable,  but  that  for  the  original  differential  equation  (1),  (4)  the  Green's  function  is 
unavailable. 
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Theorem  2.8  Suppose  $  :  [a,c]  — ►  R”  is  the  unique  solution  to  (■‘i),  (4).  Suppose  further  that 
po  :  [a,c]  — ►  L(R"^”)  is  continuous,  and  To  ;  [a,r)  — »  L(R"*‘”)  is  a  fundamental  matrix  for 
the  equation 

$'(a:)  +  Po(x)-$(i)  =  0,  (29) 

and  Gq  :  [a,c]  X  [a,c]  — *  L(R"*")  is  the  Green’s  function  for  the  boundary  value  problem  (25), 
(4).  Suppose  finally  that  the  matrix  Dq  defined  by  the  formula 

Do  =  •  To(a)  +  C  •  To(c)  (30) 

is  nonsingular.  Then  $  can  be  obtained  via  the  formula 

$(2:)=  [  Goix,t)-a(t)dt,  (31) 

J  a 

with  a  :  [a,c]  — >  R”  the  solution  to  the  second  kind  integral  equation 

oix)-poix)-J  Goix,t)-<r{t)dt  =  F^J  Go{x,t)  ■  a{t)  dt.x'^  .  (32) 

Proof.  Since  Do  is  nonsir.gular,  the  background  Green’s  function  Go  is  unique,  and  there¬ 
fore  $  can  be  obtained  from  (31).  Now,  (32)  is  obtained  by  substituting  (31)  into  (5).  □ 

2.3.  Green’s  Functions  for  Particular  Equations 

Lemmas  2.1  -  2.4  of  this  subsection  provide  fundamental  matrices  and  Green’s  functions 
for  two  particular  types  of  boundary  value  problems.  Lemmas  2.1,  2.2  are  easily  verified  by 
substituting  formulae  (34),  (35)  into  (7),  (22).  Similarly,  Lemmas  2.3,  2.4  arc  verified  by 
substituting  formulae  (37),  (38)  into  (7),  (22). 

Lemma  2.1  A  fundamental  matrix  Tq  for  the  equation 

$'(x)  =  0  (33) 

is  given  by  the  formula 

To(x)  =  /„,  (3-1) 

with  n  the  dimensionality  of  the  problem  (33),  and  x  G  [fl.c]  (in  accordance  with  standard 
practice,  /„  denotes  the  unity  operator  R"  — *  R"^. 


Lemma  2.2  The  Green’s  function  Go  corresponding  to  the  equation  (33)  subject  to  boundary 
conditions  (4)  is  given  by  the  formula 


Go(x,t) 


-  {A  F  C)-^  ■  C  (/  <x), 
-{A  +  C)-^-C  (t>x). 


(35) 
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(36) 


Lemma  2.3  For  any  X  £  R,  a  fundamental  matrix  To  for  the  equation 

$'(i)  +  A  •  $(a:)  =  0 


is  given  by  the  formula 


To(a:)  =  e-^"-/„, 


with  n  the  dimensionality  of  the  problem  (36),  and  x  £  [0, 1]. 


Lemma  2.4  The  Green’s  function  Gq  corresponding  to  the  equation  (36)  subject  to  boundanj 
conditions  (4)  is  given  by  the  formula 


.  (4  ^  ■  C)-^  -C  (i  <  x). 

•(/! +e-^ -C)-! -C  (t>x). 


(3S) 


2.4.  Linear  Transformations  for  Problems  with  Singular  Do  or  Dy 
The  purpose  of  Theorem  2.9  is  to  permit  the  conversion  of  a  problem  (1),  (3)  to  a  second 
kind  integral  equation  (28).  For  most  problems,  Theorems  2.3  and  2.7  allow  such  a  conversion, 
but  Theorem  2.3  cannot  be  used  when  the  matrix  Dfj  defined  by  (9)  is  singular,  while  Tlieorein 
2.7  cannot  be  used  when  the  matrix  Dq  defined  by  (26)  is  singular.  We  remove  these  obstacles 
in  this  subsection  by  providing  a  scheme  which  reduces  a  problem  of  the  form  (1).  (3)  witli 
singular  matrices  Dq,  Ds  to  a  problem  of  the  same  form  with  nonsingular  Dq.Ds- 

Theorem  2.10  generalizes  Theorem  2.9;  it  permits  the  conversion  of  nonlinear  problems  of 
the  form  (5).  (0)  to  nonlinear  integral  equations  of  the  form  (32). 

Remark  2.3  If  only  the  matrix  Do  is  singular,  one  can  always  choose  a  new  background 
Green’s  function  Gq  for  which  Dq  will  be  nonsingular.  However,  we  have  found  that  foi  most 
problems  it  is  easier  to  develop  a  transformation  of  the  type  described  in  this  subsection  than 
to  develop  an  alternate  background  Green’s  function. 

Theorem  2.9  Suppose  $  :  [a,c]  — *  R”  is  the  unique  solution  to  the  problem  (1).  (4)-  Suppos( 
further  that  Tq  :  [n,c]  —  L(R’”‘")  is  a  fundamental  matrix  for  the  background  equation  (4o). 
Suppose  finally  that  there  exists  :  [n,c]  —  L(R”’‘”)  such  that  'I'  £  C’[«,c],  dot  'l'(.i  )  ^  0  for 
all  X  £  [n,c],  and  the  matrix 

Do  =  A-^{a)-roia}  +  C  -^(cj-Tolc)  (39) 

is  nonsinrjuinr.  Then  the  equation 

r'(a-)  +  'i'-^x)  •(^''(T)  +  p(r)-  'i'(x))-r(x)  =  ^'-^(x)  •  /(x).  (-10) 

subject  to  boundary  conditions 

A-<l!(a)-T(a]-\-C  ■'<l/{c)-T{c)  =  0.  (-11) 

has  a  unique  solution  F  :  [a,c]  — ►  R”,  and 

$(x)  =  ^'(x)-r(x),  (42) 

for  all  x  £  [a,  c]. 
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Proof.  We  immediately  obtain  (41)  by  substituting  (42)  into  (4).  Now,  substituting  (12) 
and  its  derivative  into  (1),  we  get 

'^'(i)  ■  r(x)  +  ^{x)  ■  r'(x)  +  p(i)  .  5»(ar)  •  r(x)  =  fix),  (43) 

and  obtain  (40)  by  combining  (43)  with  the  fact  that  '^(x)  is  nonsingular  for  all  x  G  [a.c].  □ 

Remark  2.4  Clearly,  the  transformed  problem  (40),  (41)  satisfies  the  conditions  of  Theorem 
2.7,  as  Do  defined  by  (39)  is  nonsingular.  However,  for  many  problems,  the  boundary  condition 
matrix  defined  by  the  formula 

Df,  =  A  ■<^(a)-¥C -^(c).  (44) 

is  singular,  and  therefore  the  transformed  problem  fails  to  satisfy  the  conditions  of  Theorem  2.3. 
If  one  needs  to  use  the  results  of  both  The  rem  2.7  and  2.3,  one  must  choose  a  transformation 
such  that  both  Dq  and  are  nonsingular. 

Of  course,  it  is  easier  to  choose  transformations  when  Do  =  This  is  true  when 

the  background  Green’s  function  is  chosen  to  correspond  to  the  equation  =  0.  By  Lemma 
2.1,  the  fundamental  matrix  for  this  equation  is  T  =  /„;  the  equivalence  for  this  fundamental 
matrix  of  (39)  and  (44)  is  readily  apparent. 

Theorem  2.10  is  the  nonlinear  analogue  of  Theorem  2.9;  the  proofs  of  the  two  theorems  are 
nearly  idejitical. 

Theorem  2.10  Suppose  $  :  (fl,c]  —  R'‘  is  the  unique  solution  to  the  problem  (5).  (4)-  Suppose 
further  that  Tq  :  [a,c]  — ►  Z,(R"^”)  is  a  fundamental  matrix  for  the  background  equation  (25). 
Suppejsc  finally  (hat  there  exists  :  [a,c]  — *  LfR"**”)  such  that  ’J'  £  C^[a,c],  det  i>{2')  ^  0  for 
all  X  £  [a.c],  and  the  matrix 

Do  =  A-  'l'(a)  •  To(«)  +  C  •  'I'(c)  To(c)  (4-3) 

is  nonsingular.  Then  the  equation 

r'(x )  +  (x  )'I''(x)  •  r(i)  =  ^'“'(x)  •/'(  4'(x)  •  r(.r),  .r)  (4G) 

subject  to  boundary  conditions 

4  •  ^*(0)  ■  r(a)  +  C  •  \l'(c)  •  r(c)  =  0.  (47) 

has  a  unique  solution  F  :  [a,c]  — >  R",  and 

4>(x)  =  ^'(x)  •  r(x),  (48) 


for  all  X  £  [a,  c]. 
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2.5.  Newton’s  Method  for  Nonlinear  Boundary  Value  Problems 
Theorems  2.4,  2.8  of  Section  2.2  reduce  nonlinear  boundary  value  problems  of  the  form 
(5),  (6)  to  nonlinear  second  kind  integral  equations  of  the  form  (32).  In  this  subsection,  \vc 
describe  the  convergence  properties  of  the  well-known  Newton’s  method  as  applied  to  the 
latter  (Theorem  2.12),  and  reduce  each  step  of  Newton’s  algorithm  to  the  solution  of  a  linear 
boundary  value  problem  of  the  form  (1),  (4)  (Theorem  2.11). 

Theorem  2.11  permits  each  Newton  iterate  Sk  defined  by  (11)  to  be  expressed  as  the  solution 
to  a  second  kind  integral  equation. 


Theorem  2.11  Suppose  K  :  — *  R"  in  C’'[a,c]  defined  by  (10)  is  Frechet  differentiable 

at  every  point  x  £  [a,c],  and  4>/;  :  [u,c]  — >  R"  is  defined  for  all  k  =  O.l. .. .  via  the  formula 

$^.(T)=  r  Go{xJ)-ak(1)d1.  (49) 

J  a 

with  Ok  :  [n.c]  — >  R”  defined  by  (12),  and  Cq  :  [u,c]  X  [fi,cj  — •  L(R’”‘")  the  Green's  function 
fur  (25).  (4).  Then  the  Newton  iterates  :  [fl.c]  — -  R'*  €  C^fn.cj  given  by  Definition  2.1) 
satisfy  the  equation 

h{i]  +  Dk{x]-  r  Go(xJ).6k{t)dt  =  gDx)  (oO; 

J  a 

for  a!!  k  =  0,  1 , .  .  .,  with  ilk  :  [a,  c]  —  L(R"’‘’’ )  defined  for  alt  k  =  0.  \ by  the  furnnilei 


Qkix) 


dF{^k{x),x) 

d^k 


Po(x). 


(ol) 


and  Qk  :  [a,c]  — >  L(R’^’‘")  defined  for  all  k  0,1, ...  by  the  formula 


gk{x)  =  po(x)  ■  ^k(x)  +  F(^Ax),x)  -  ak{x). 


Proof.  (50)  is  obtained  by  substituting  the  Frechet  derivative  of  the  function  J\  into  (12). 
and  substituting  (49),  (51),  (52)  into  the  resulting  equation.  □ 

The  coiv.'ergcnce  properties  of  .Newton’s  metliod  have  been  thoroughly  studied.  Theorem 
2.12  is  one  fundamental  result,  and  can  be  found,  for  example,  in  [12]  (in  a  slightly  dificicnt 
form ). 


Theorem  2.12  Supjmsc.  is  the  unique  solution  to  (5).  ()),  b  is  the  solution  to  (22).  and  k 
IS  the  unique  solution  to  (50)  (so  that  the  linearization  (50)  to  the  equation  {32}  is  nonsingulav 
at  a).  Then  there  exists  c  >  0  such  that  for  any  (Tq  :  [a,c]  —  R”  satisfying  the  condition 

Iko  -  ^11  <  < 

and  Newton  iterates  Ok  :  [u,c]  —  R"  defined  by  (11), 

1.  \\ak  —  d||  <  f  for  all  k  =  1,2,..., 

2.  lim  Ok  ~  b, 
fc— *00 

3.  Ok  conicrges  to  b  quadratically. 
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2.6.  A  Lemma  from  Linear  Algebra 

Given  a  perturbation  of  the  unity  operator  '■  the  Lemma  2.5  provides 

its  inverse.  It  is  normally  used  when  the  rank  of  the  perturbation  is  low,  is  a  particular  case  of 
the  Sherman-Morrison  formula  (see,  for  example,  [9]),  and  is  easy  to  verify  directly. 

Lemma  2.5  For  any  two  vectors  U,V  6  such  that  U  ^  In- 

^  +  U  (In-  U)-^  ■  v^.  (54) 

III.  The  Analytical  Apparatus 

In  the  remainder  of  this  paper,  we  assume  that  the  solution  to  the  problem  (1),  (4)  is 
being  sought  on  the  interval  [a,c],  and  that  6  is  some  intermediate  point  (a  <  b  <  c).  The 
fundamental  observation  on  which  Algorithm  A  is  based  is  that  the  solution  to  the  integial 
equation  (2S)  on  the  entire  domain  [a,c]  can  easily  be  constructed  from  the  solutions  of  two 
independent  integral  equations,  one  defined  on  [a,  6]  and  one  on  [i».c].  This  leads  naturally  to 
a  recursive  algorithm,  in  which  independent  solutions  on  a  large  number  of  subintcrvals  are 
successively  merged  until  the  fuU  solution  is  obtained.  A  precise  formulation  of  the  construction 
and  the  resulting  numerical  scheme  will  require  some  notation. 

3.1.  Notation  We  will  denote  the  subintervals  [o.6]  and  [6,c]  of  [fl,c]  by  A  and  B.  respec¬ 
tively.  For  convenience,  we  write  the  integral  equation  (28)  in  the  form 

cr{x)  +  p(x)-  [  Go(Xyt)  ■  o(t)dt  =  fix).  (55) 

J  a 

with  fix)  —  p(j)  -  poix).  and  Gq  :  [a,c]  x  [a,c]  —  L(R"’‘")  the  background  corresponding  tn 
the  equation  (25)  subject  to  boundary  conditions  (4). 

We  define  the  operator  P  :  (X^[o,c])"  —  (X^[a,c])”  corresponding  to  (55)  by  the  formula 

P((7)(a  )  =  a(j)  +  p(x)  ■ 

so  that  we  have 

Po  =  f. 

M’e  will  renjuim  the  four  operators 

Paa  ■■  (X2[a,6])"^(X2[„,i])"  , 

Pab  ■■  (L"[6.c])"-(X^[n,fi])"  , 

Pba  :  (X'[a,6])'‘^(X2[fi,e])"  , 

Pbb  :  (X'[h,c])"-(X'(6,c])", 


r 


Go(x.t)  ■  a[1)(ll. 


'  oCi  I 
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defined  by  the  formulae 


Paa{<^){x) 

Pab{^){x) 

Pba{<^){x) 

PBB{<y){x) 


fl> 

a{x)-\-p{x)-  j  Go{x,i)- a{t)dt, 
Ja 

p{x)-J  Go(x,t)-a(t)dt, 

fb 

p{x)-  /  Goix,t)-  cr{t)dt, 

Ja 

a{x)  +  pix)-j  Go{x,t)  ■  a{t)di. 


We  define  the  operator  Q  :  iL^[a,c])''‘^^  —  (l2[a,c])”’‘"  by  the  expression 

Q{x)(x)  =  xix)  +  p(x)-  I  Goix,t)-  x{t)dt. 

Ja 


We  additionally  require  the  four  operators 


Qaa 

Qab 

Qba 

Qbb 


(50) 


(GO) 

(61) 


(02) 


defined  by  the  formulae 

Q^,4(\)(.r) 

QABix)ix) 

Qba(x){x) 

QBB(\)ix) 


Xix)  +  p{x)-  t  Go{x,1)-x(ndt, 

J  a 

p{x)-  Go{x,t)-x(t)dt. 

Pix}-  f  Go(xJ)-x{t)dt, 

Ja 

Xix)  +  P(x)-J^  Goix,i)  ■  \U)  dt. 


(03) 


(64) 


(65) 


(66) 


We  also  require  the  functions  L(R'”‘")  defined  by  the  formulae 

i’(x)  =  p(i)-To(x),  (60 

t'JO  =  To'(0  +  -/o(0, 

v„{t)  =  ^(0, 

with  To  the  fundamental  matrix  for  equation  (25),  and  Jq  :  [a,c]  —  L(R”’'")  defined  liy  the 
formula 

^(0  = -^o’ •C-To(c)-To’(0,  ('6) 

with  the  matrix  Dq  defined  by  (26),  and  the  matrix  C  given  by  (4). 
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Given  a  function  /  e  (L^[a,c])”,  we  will  follow  the  convention  of  denoting  its  restriction 
to  A  and  B  by  /|^  and  /jg,  respectively.  Similarly,  given  a  function  0  G  we  will 

denote  its  restriction  to  A  and  B  by  4>\a  ^nd  V'iSi  respectively.  Assuming  that  the  operators 
PaAiPbb  are  nonsingular,  we  define  the  functions  t]a  :  A  R”,  Tjs  ■  B  R”  via  the  formulae 


TlA 

=  PA\if\A). 

(71) 

VB 

=  -PbbC/ib). 

(72) 

Similarly,  assuming  that  the  operators  QiQaAiQbb  are  nonsingular,  we  then  define  the 

map- 

pings 

X 

:  [a 

<f>A 

:  A 

^  L(R"'‘”), 

4>b 

:  B 

^  LCR”*"”), 

via  the  formulae 

X 

= 

(73) 

4>a 

=  QA\i^\Ah 

(7-1) 

<i>B 

~  Qbb(^\b)- 

(75) 

Finally,  we  will  define  six  matrices  o^,q^,  of  ,af,  €  L(R"*‘")  by  the  formulae 

= 

1  VL{t)  ■  4>A{t)  dt  , 

/a 

(76) 

= 

1  Vii{t)-4>A{t)dt  , 

J  a 

(77) 

= 

^  vdt)  ■  <l>Bit)  d1  , 

(7S) 

= 

^  Vnit)  ■  <l>Bit)  dt  , 

(79) 

Oil 

= 

/  ■  xit)  di  , 

Ja 

(SO) 

= 

f  v,i(t)  ■  x{t)  dt  , 

Ja 

(SI) 

and  six  vectors 

e  R” 

via  the  formulae 

= 

1  Vi,{t)-r)A{t)dt, 

Ja 

(S2) 

= 

f  Vn(t)-'qA(1)dt, 

Ja 

(S3) 

= 

VL{i)-m{t)dt, 

(84) 
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= 

VKit)-T]Bii)dt, 

(85) 

h  = 

f  vj,{t)-a{i)dt, 

Ja 

(80) 

6fi  = 

f  VRit)-(T(t)dt, 

Ja 

(87) 

with  a  the  solution  to  equation  (57). 

3.2.  Analysis  of  the  operators  Pab,Pba 

In  this  subsection,  we  observe  that  each  of  the  operators  Pab  Pba  is  of  rank  n,  and 
give  simple  expressions  for  these  operators. 

Lemma  3.1  In  the  notation  of  the  preceding  subsection, 

Pab  =  i’lA  ■ 

Pba  =  i’\B-^I- 

Proof.  We  obtain  (88)  by  observing  that  x  <  t  for  any  x  €  E  [i'.c],  and  by  using 

(59)  and  (69).  Similarly,  (89)  follows  from  the  combination  of  (60),  (69)  and  observing  that 
X  >  t  for  any  x  6  [h,  c],  i  €  [a,  b].  ° 

3.3.  Recursive  solution  of  the  integral  equation  (57) 

We  now  consider  the  original  integral  equation  (57) 

Pa  =  /. 

The  main  result  of  this  subsection  is  the  following  lemma,  whicli  constructs  tlie  soliitioii  o 
of  equation  (57)  from  VAiVB  of  equations  (71)  and  (72). 

Lemma  3.2  If,  in  the  notation  of  Section  3.1,  all  six  operators  P,  Paa,  Pbb  Q-  Qaa-  Qbb 


are  nonsingular,  and  the  matrices 

Aj ,  Aj  G  L(R" 

‘’‘")  defined  by  the  formulae 

0 

1 

11 

< 

a®. 

( 90 

A 2  —  In  —  Or  ' 

(91 

art  also  nonsingular,  then 

(^\A 

T]A  +  <f>A-  AJ*  • 

(«f  ■ 

(92 

= 

VB-f  <Pb  ■  AJ"’ 

(93 
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Proof.  Using  definitions  (56)  -  (61),  the  integral  equation 

P<T  =  f 


caji  be  rewritten  in  the  form 


Paa{(^\a)  +  PABio^B)  =  f\A-, 

Pba{o\a)  +  Pbb(CT|b)  =  /|B-  (^5) 

The  expansions  (88)  and  (89)  for  Pab  a-nd  Pba-,  respectively,  can  then  be  used  to  obtain  an 
explicit  solution  to  the  coupled  equations  (94)  and  (95)  in  terms  of  the  functions  t]a,  rjs,  4>a- 
and  (pB  defined  by  (71),  (72),  (74),  (75),  respectively.  Indeed,  applying  the  operator  P^\  to 


equation  (94)  and  the  operator  Pgg  to  equation  (95),  we  have 

^\A  +  PaA  ■  PAB{(r\B)  =  PAAif]A)^  (9G) 

PbB  '  PsAi(^\A) ^\B  =  PshiflB)-  (^”) 

Substituting  (88)  and  (89)  into  (96)  and  (97)  yields  the  formulae 

(^\A  + Paa-‘^\a-'"1-(^\B  =  Va,  (9^) 

Pbb  ■  '^'\B  ■  vj  ■  o\A  +  o^B  =  7?b,  (99) 

or 

<r\A  +  <t>A  ■  ■  o^B  =  (100) 

<Pb -vl o^B  =  VB,  (101) 

where  we  have  used  the  definitions  (74),  (75)  for  4>a  and  4>Bi  respectively.  Now,  multiplying 
(101)  by  4>a  ■  and  subtracting  it  from  (100),  we  obtain 

(/(B2)n  -  <Pa  ■  ■  <f>B  ■  '(’I)-<^\A  =  Va  -  <Pa  ■  ■  rjB-  (102) 

Similarly,  multiplying  (100)  by  4>b  •  and  subtracting  it  from  (101)  results  in  the  equation 

(I^L2)n  -  4>b  ■  vl  ■  4>a  ■  v^)  ■  <^\B  =  riB  -  <Pb  ■  vl  ■  lA-  (103) 

Due  to  (76),  (79),  and  (82),  (85)  we  can  rewrite  these  equations  in  the  form 

(/(B2)n  -  4>a  ■  (of  •  vI))-<7\a  =  Va-  <t>A-f>R-.  (104) 

{I(l2)n  -  <Pb  ■  (a^  •  v'^))■a^B  =  VB  -  <Pb  ■  (105) 

By  application  of  Lemma  2.5,  we  obtain 

(^\A  =  +<Pa  -ih  -  Or  -vJ  -(pA^^  -  of  ■vl)-(TtA  -  <pA  '^r),  (lOG) 

OlB  =  +<f>B  -  {In  -  Oi  •  'oI-<PbT^  •  vl)  ■  {t)b  -  <Pb  )•  (107) 

The  equations  (92),  (93)  are  now  obtained  from  equations  (106),  (107)  and  equations  (90). 
(91).  □ 


Remark  3.1  Suppose  that  bi  and  62  are  a  pair  of  real  numbers  such  that  a  <  61  <  62  <  c,  and 
that  the  interval  [hj,  62]  is  denoted  by  C.  We  will  denote  by  Pcc  the  restriction  of  the  operator 
P  to  the  interval  C,  and  denote  by  Qcc  the  restriction  of  the  operator  Q  to  the  interval  C. 
Assuming  that  Pcc, Qcc  are  nonsingular,  we  define  the  functions  rjc  '■  C  -*  R",(5!>c  :  C  — 
L(R">'")  by 

VC  =  peMAc),  (lOS) 

4>c  ~  Qcci^\c)  •  (109) 

By  applying  the  above  lemma  twice  (once  for  the  subinterval  [a,ii]  and  once  for  [a,  62])! 
may  easily  observe  that  there  exists  A  6  R"  such  that 

(t(x}  =  T}c{x)  +  <f>c(x)  ■  X  (110) 

for  all  X  €  C.  The  exact  expression  for  the  vector  A  is  complicated,  but  irrel^^vanf  for  the 

purposes  of  this  paper.  The  existence  of  a  relation  of  the  form  (110),  however,  will  be  critically 

important  in  Section  4. 

3.4.  Further  Analytical  Results 

We  now  collect  a  number  of  identities  which  are  necessary  for  Algorithm  A,  to  be  prescjited 
in  Section  4.  Corollary  3.1  is  similar  to  Lemma  3.2,  but  uses  the  matrix  valued  function  V' 
in  place  of  the  vector  valued  function  /  to  obtain  an  analytical  expression  for  the  function  \ 
defined  by  (73). 

Corollary  3.1  If,  in  the  notation  of  Section  3.1,  all  six  operators  P,  Paa,  Pbb,Q,Q  .\a,Q  BB 
are  nonsingiilar,  then 

X|4  =  <t)A-  (111) 

X|B  =  </-B-Ar'.(/„-Q^).  (112) 

with  the  matrices  and  of  defined  by  equations  (76)  and  (79),  and  the  matrices  Ai.Aj 
defined  by  equations  (90)  and  (91). 

Proof.  Substituting  in  equations  (106),  (107)  the  functions  (()a.4>b  defined  by  (71).  (7')) 
for  the  functions  t]a,Vb  defined  by  (71),  (72),  and  the  matrices  Q^,af  defined  by  (76).  (79) 
for  the  vectors  1^,6^  defined  by  (82),  (85),  we  obtain 

Am  =  <t>A  -  4>a  ■  +  4>a  ■  ■  of  ■  a'^  -  (t>A  ■  ■  af  ■  ,  (113) 

X|B  =  4>B  -  <^B -O^  +  epB  ■  ■o'^  -of  ~  4>B  ■  -O^  -of  ■  q'^.  (114) 

The  expressions  (111),  (112)  are  now  easily  obtained  from  the  equations  (113),  (114).  □ 

We  w'ill  also  require  analytical  expressions  for  the  inner  products  6^  and  Sr  defined  by  (86). 
(87)  in  terms  of  the  restricted  inner  products  Sf,Sf,6^  and  6^  defined  by  (82)-(85). 
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Corollary  3.2  If,  in  the  notation  of  Section  3.1,  all  six  operators  P,Paa,  Pbb,QiQaa\Qbb 
are  nonsingular,  then 


h  =  j  Vi{t)-(T{t)dt  =  j^V],{t)-a^A{l)dt-\- Vj,{i)-0)B{l)dt 

=  +  Sf  +  at  •  “(of  ■  6t  -  «l)  + 


(115) 


6ft  -  j  Vft{t)-a{t)dt  =  j^Vn{t)-cr\A(t)dt-\- Vn{t)-a\^B{l)di 

=  ^ J  +  a^  Aji  “(a®  .  6t  -  )  +  of  •  Ar'  •  (of  •  -  6,^  )•  ( 116) 

Proof.  Multiplying  equation  (92)  by  vj  and  vj,  and  equation  (93)  by  vj  and  uj,  we  obtain 


Ja 

^  i’i(0 

*/  a 


+  -Aj^  -(of -St  -6^), 
+  of  .A^'•(a^6f -6,^), 

,5^  +  o^  •  Aj^  -(af  -<5^  -  ), 

+  «f -Ar'  -^t). 


(117) 

(US) 

(119) 

(120) 


Now,  expressions  (115),  (116)  are  easily  obtained  from  (117)-(120).  63 

CoroDary  3.3  is  similar  to  Corollary  3.2,  but  uses  x,  the  matrix  valued  function  defined 
by  (73),  in  place  of  a,  the  vector  valued  function  defined  by  (57).  While  the  two  corollaries 
concern  different  objects  (the  vectors  in  Corollary  3.2,  the  matrices  0^,0^  in  Corollary 

3.3),  their  proofs  are  nearly  identical. 


Corollary  3.3  If  in  the  notation  of  Section  3.1,  all  six  operators  P,  Paa,  Pbb^Q  aa-Q  BB 

are  nonsingular,  then 


Ql 


[  Vf,{t)-x{t] 
Ja 

*^L  ■  A2  '  (dn 


dt  =  J  VLit)  ■  X\Ail)  dt  +  Vfit)  ■  X\Bi^)dt 
~  )  "b  ■  Aj  ‘  {dn  ~  ^l)i 


(121) 


Or  =  /  Vh{t)-X{l)dt-  !  Vft{t)-x\Ail)dt+  I  Vn{t)  ■  X\Bil)  dt 

Ja  Ja  -tfc 

=  Or  •  AJ*  •  (/n  -  of  )  +  Or  •  Aj  ’  •  (/n  -  of  )• 

Finally,  combining  Lemma  3.2  with  the  expressions  (111)-(112),  we  have 


(122) 
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Corollary  3.4  Suppose  that  in  the  notation  of  Section  3.1,  all  six  operators  P,  Paa<  Pbb- 
QiQaAjQbb  O’"®  nonsingular.  Suppose  further  that  the  function  F  is  defined  by  the  formvlu 


-F’(i)  =  X  •  (123) 

with  X  6  R".  Then  on  the  interval  [a,i], 

F{x)  =  d>A(x)- n  + t]a(x),  (12.I) 

with  p,  £  R"  defined  by  the  formula 

p  =  A^\X-af.iX-6t)-6^).  (123) 

Similarly,  on  the  interval  [6,  c], 

F(x)  =  cI)b{x)  ■  1/ +  t)b(x),  (12G) 

with  u  £  R"  defined  via  the  formula 

o  =  ■{X-ai-(X-6f)-6-^).  (127) 

Proof.  Restricting  (123)  on  the  subintervals  A,B  of  [fl,c],  respectively,  we  have 

^\A  =  +  (12S) 

f\B  =  XlB'-^  +  ^^lB-  (129) 

Combining  (12S),  (129)  with  (92),  (93),  (111),  (112),  we  obtain 

F\a  =  (t>A-^2^-iJn-of)-X  +  iriA  +  d>A-^2^-iaf-6-^~S^)),  (130) 

Fjs  =  d>B  ■  ■  (In  -  Q'l)- ^  + (jlB  +  d>B  ■ -{at  ■  -  Fl))-  (1-11) 

Now.  the  expressions  (125),  (127)  immediately  follow  from  the  comparisons  of  (12-1).  (12G)  with 
(130).  (131 ),  respectively.  □ 


IV.  Description  of  the  Algorithms 

We  turn  now  to  the  construction  of  the  fast  algorithm  for  the  solution  of  the  inlogial 
equation  (57) 

Pcr  =  f  , 

based  on  the  apparatus  developed  in  Section  111.  The  main  tool  at  our  disposal  is  the  aljility 
to  merge  the  solutions  of  restricted  versions  of  the  integral  equation  in  adjacent  subintervals 
(Lemma  3.2).  As  this  suggests  a  recursive  procedure,  we  begin  by  subdividing  the  wliole 
interval  [a,c],  on  which  the  solution  to  (57)  is  sought,  into  a  large  number  of  subintervals.  For 
the  sake  of  simplicity,  we  assume  that  m  is  a  positive  integer  and  that  M  =  2’”  is  the  nunibei 
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of  subintervals  created.  The  boundary  points  of  the  subintervals  are  then  defined  by  a  strictly 
increasing  sequence  of  numbers 

(132) 

with  bi  =  a  and  6a/+i  =  c.  For  each  t  =  1, . . Af,  we  define  the  interval  via  the  expression 

5r  =  [6.,i.+i],  (133) 

and  create  a  hierarchy  of  intervals  by  recursively  merging  adjacent  pairs.  That  is,  for  eacli 
j  =  m  —  1, . . . ,  1,0,  and  i  =  1, . . .,  M,  we  define 

=  (134) 

We  will  refer  to  each  fixed  /  as  a  level.  We  will  also  refer  to  the  two  intervals  B2t-i  -^2^' 

as  children  and  to  the  larger  interval  fl-  as  a  parent. 

It  is  obvious  that 

Bi  =  (135) 

and  that  for  each  level  /, 

2* 

[fl,c]=ijB'.  (136) 

t=i 

4.1.  Notation 

Generalizing  the  notation  of  Section  III,  we  will  denote  by  Pij  the  restriction  to  the  interval 
5-  of  the  integral  operator  P,  so  that 

P,j(a){x)  =  a{x)  +  p{x)-  f  Go{x,t)  ■  a{t)dt  (137) 

a”*-' 

for  any  a  e  L^{B\Y.  Similarly,  we  will  denote  by  Qii  the  restriction  to  the  interval  B\  of  the 
integral  operator  Q,  so  that 

Qi,l{x){x)  =  xi^)  + Pi^)-  I Go{xJ)-xU)dt  (13S) 

■'^l+(t-l)  2”'-‘ 

for  any  x  €  L'^{B\Y'^'^.  For  each  B\  we  will  define  the  functions  t;,,/  :  B[  —  R",p,,;  :  B[  — 
L(R"''")  as  the  solutions  of  the  equations 

P.,;(^,,/)  =  /|fi.,  (139) 

Qi,l{4‘i,l)  — 

provided  the  operators  Pi,i,Qi,i  are  nonsingular. 


(140) 


Remark  4.1  Suppose  now  that  the  operators  are  nonsingular  on  the  interval  D[. 

Then,  due  to  (110),  there  exists  X''‘  6  R"  such  that 


<7(1)  = 

(141) 

for  all  I  G  BI- 

For  each  /  =  0, 1, . . 

.,m,  and  i  =  1,2,..., 2*,  we  define  the  matrices  a'A , 

0/  G  L(R"’'")  by 

the  formulae 

■  4>,j{t)dt, 

J'"-' 

(142) 

(It, 

(143) 

and  the  vectors 

G  R’^  by  the  formulae 

S'/  = 

(144) 

s/  =  ■r),j{t)dt. 

2">-' 

(145) 

4.2.  Discretization  of  the  Restricted  Integral  Equations 
Choosing  an  integer  p  >  1,  we  construct  the  p  scaled  Chebyshev  nodes 


cos 


(2j  -  l)7r 


2p 


+ 


j  =  1,2 . p 


(MG) 


on  each  of  the  intervals  5'",?  =  1,2,...,M.  We  then  discretize  the  two  integral  equations  ( 139). 
(140)  via  a  Nystrom  algorithm  based  on  p-point  Chebyshev  quadrature  (.'^ec.  for  example.  [11])- 
The  resulting  approximations  to  the  functions  r],j,  4>,j  at  the  nodes  t-  will  be  denoted  by 


fl>,t  =  , 


respectively. 

Remark  4.2  It  is  well-known  that  the  order  of  convergence  of  the  approximations  fhj-Oi.i  to 
the  functions  is  p.  Since  all  subsequent  steps  in  the  construction  of  an  approximate 

solution  a  to  the  integral  equation  (57)  are  analytic,  the  convergence  rate  of  the  full  algorithm 
depends  entirely  on  the  parameter  p.  For  example,  by  using  16  scaled  Chebyshev  points  on 
each  subinterval  at  the  finest  level,  one  obtains  a  sixteenth  order  method. 
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Remark  4.3  The  algorithm  of  this  section  makes  extensive  use  of  the  apparatus  of  Chebyshev 
interpolation,  quadratures,  composite  quadratures,  etc.  This  apparatus  is  quite  well-developed, 
and  can  be  found  in  various  forms  in  [6],  [8],  [10].  For  a  detailed  description  in  the  form  most 
convenient  for  our  purposes,  we  refer  the  reader  to  [11], 


4.3.  Informal  Description  of  the  Algorithm  for  Linear  DDEs 

We  begin  by  directly  solving  the  two  integral  equations  (139),  (140)  on  each  subinterval  5"’ 
at  the  finest  level,  as  discussed  in  the  preceding  subsection.  Equation  (141)  then  shows  that 
a  restricted  to  BJ"  can  be  expressed  as  a  linear  combination  of  the  two  solutions  77, 

Thus,  it  remains  only  to  determine  the  coefficients  A‘’"‘  6  R”  for  each  of  the  M  subintervals 
BY'.  Fortunately,  this  can  be  done  recursively.  To  see  this,  suppose  that,  at  some  coarse  level 
/  <  m  -  1,  we  are  given  the  coefficient  for  the  subinterval  B[.  Then  Corollary  3.4  provides 
formulae  for  the  calculation  of  the  corresponding  coefficients  £  R"  for  the  two 

child  intervals  Bjlii  and  respectively.  On  the  coarsest  level,  we  observe  that  A“  *  =  0. 

i.e.  the  solution  of  equation  (139)  on  the  whole  interval  [a,c]  is  simply  a. 

However,  the  formulae  (125)  and  (127)  of  Corollary  3.4  contain  the  matrices 

and  the  vectors  These  quantities 

are  also  computed  recursively  but  in  the  opposite  direction,  namely,  from  the  finest  level  to  the 
coarsest.  They  are  certainly  available  at  level  m  directly  from  the  definitions  (142)-(145).  For 
the  interval  B;  at  any  coarser  level  I  <  m  —  1,  Corollaries  3.2  and  3.3  describe  liow  matrice.s 
and  vectors  are  obtained  from  the  matrices  ai,Qfi  and  vectors  61,6^  of  the  two 

child  intervals. 

To  summarize,  the  algorithm  consists  of  three  parts.  First,  a  sufficiently  fine  subdivision 
61, 62,  •  •  of  the  interval  [a,c]  is  chosen  so  that,  on  each  of  the  intervals  Bj,mi  the  functions 

can  be  accurately  represented  by  a  low  order  Chebyshev  expansion.  On  each  of  the 
intervals  B,,„;,  the  equations  (139)  -  (140)  are  solved  (approximately)  by  direct  inversion  of 
the  linear  system  arising  from  a  Nystrom  discretization.  Second,  the  matrices  and 

vectors  6'i,6'k  are  computed  in  an  upward  sweep,  beginning  at  the  finest  level  m.  Finally, 
the  coefficients  A'’^  are  computed  in  a  downward  sweep,  beginning  at  the  coarsest  level.  The 
desired  function  o  is  then  recovered  on  each  subinterval  from  equation  (141). 

The  following  is  a  more  detailed  description  of  the  numerical  procedure. 

Algorithm  A 


Comment  [Define  the  computational  grid.) 

Create  M  =  2"’  subintervals  on  [a,  c]  by  choosing  a  sequence  of  boundary  points  61 , 62.  •  ■  ■  >  ,  ^AZ  +  i 

with  6i  =  a  and  6a/4i  =  c.  Choose  the  number  p  of  Chebyshev  nodes  on  each  interval  B,"*  =  [6,,  6,  +  i] 
for  i  =  I, . . . ,  M .  Determine  the  locations  of  the  scaled  Chebyshev  nodes  ,  r?, . . . ,  rf  on  each  interval 
BY',  and  evaluate  the  functions  /,  V'  at  these  nodes,  obtaining 

Step  1. 

Comment  [Construct  the  approximate  solutions  on  each  interval  BY'  ] 
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do  i  =  1,2, 

(1)  Construct  the  two  p  n  x.  p  n  linear  systems  on  B,'  obtained  through  a  Nystrom 
discretization  of  the  corresponding  integral  equation. 

(2)  Solve  the  two  p  ■  n  x  p  ■  n  linear  systems  on  B[  by  Gaussian  elimination, 
obtaining  the  values 

end  do 


Step  2. 

Comment  [Construct  the  matrices  and  vectors  on  each  interval  B,"'  at  ihr  finest 

level.] 

do  i  =  1,2,...,  A/ 

Evaluate  the  matrices  and  vectors  using  the  p-point 

Chebyshev  quadrature  formula, 
end  do 


Step  3  (Upward  Sweep). 

Comment  [Construct  the  matrices  q1'.q‘w'  and  vectors  for  all  interval.*^  at  all  coar.s-r  l-v.  U 

/  =  m  -  l.rn  -  2 . 0.] 


do  1=  in-l,  0.  -1 
do  i  =  l.  2* 

Compute  the  matrices  a],' 
intervals 
of  Corollaries  3.2  and  3.3 
end  do 
end  do 


and  vectors 


from  the  corresponding  data  in  the  two  clnli 

■  •  21-1,/+]  x2'.'+l  1  ■  .1  I 

'n  .ot  .0;,  ),  using  the  rostili,- 


.0i  ,0^ 


Step  4  (Downward  Sweep). 

Comment  [Construct  the  coefficient  X'-'”  for  all  intervals  at  the  finest  level  ] 

Set  A'’  *  =  0. 

do  l=0,m-l 
do  i=l,  2’ 

Use  Corollary  3.4  to  compute  the  coefficients  A'"'’*  * ,  A^'  *'*'’ , 

for  the  cl’ild  intervals  B^l’  and  B^li,  from  the  coefficient  A"  '  of  the  parent  intei\al  B' 
end  do 
end  do 


Step  5. 


Comment  [Compute  the  solution  a  of  equation  (57)  at  the  nodes  r/ ,  r-, . .  . ,  rf  for  each  interval  B"‘ 
at  the  finest  level.]  '  ‘ 
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do  i=l,  M 
doj=l,p 

Determine  the  values  of  the  solution  a  of  equation  (57)  at  the  node  t\  via  formula  (141). 
end  do 
end  do 


Step  6. 

Comment  [Compute  the  solution  ^  of  equation  (1)  from  the  values  of  a.] 

Evaluate  the  integral  (27),  by  using  composite  Chebyshev  quadrature 
(see  Remark  4.5  below). 

Remark  4.4  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work  required  is  of 
the  order  ■  n^).  Step  1  involves  solving  two  (p  x  n)  x  (p  X  n)  linear  systems  for  each  of 

the  M  intervals.  Steps  2  -  5  require  no  more  than  0(Af  ■  p  ■  ■  (logp  +  n))  operations.  Since 

.V  =  M  -  p  is  the  total  number  of  nodes  in  the  discretization  of  the  interval  [a.c],  we  can  write 
the  CPU  time  estimate  in  the  forni  0(.V  -p^-n^).  The  cost  of  evaluating  the  solution  ^  of  the 
differential  equation  (55)  from  the  integral  representation  (27)  is  0[N  -logp  •  n)  (see  Remark 
4.5  below). 


Remark  4.5  The  final  s‘ jp  in  the  algorithm  involves  the  evaluation  of  an  integral  of  the  form 
(27)  at  each  of  the  Chebyshev  nodes  t-  on  each  subintc’‘val  .  namely 


/  Go{Tl,i)-a{t)dt. 
J  a 


(147) 


If  these  integrals  were  calculated  independently  for  each  r-,  the  amount  of  work  required 
would  be  of  the  form  0{N^  ■  n),  and  would  dominate  the  construction  of  the  function  In 
fact,  this  is  unnecessary,  for  we  may  write 


$(rJ)  =  T(r/; 


/  I'i-fO  •  o(t)(lt  +  f  ' 
Ja  Jb, 


Vl(1)  ■  C{1)(lt 


+  (  v„{()  ■  a{t)dt  +  f  Vf,{()  ■  cr 

■'rj  A.+  , 


{i)dt 


(  M>i 


where  we  have  used  the  representation  (22)  and  the  fact  that  lic.s  in  the  interval  B^''  = 
[f'Mf'i+i]-  Step  G  can  then  be  written  in  detail  as  follows; 


Step  6  (a). 

Comment  [Precompute  the  integrals  of  i  j.  r  and  v„  ff  on  each  subinterval  B,'"  by  Chebyshev  quadra¬ 
ture.  These  integrals  will  be  denoted  /j,  and  /„,  respectively.] 
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do  i=l,  M 

hiBn  =  f':^'vut)cTit)dt. 

end  do 


Step  6  (b). 

Comment  [March  across  interval  from  a  to  c,  computing  "I*  at  each  node  in  the  discretization.  Tlie 
variables  jTt  and  Jr,  will  be  used  to  accumulate  the  integrals  •  <T{i)di  and  Vri(t)  ■  cr(i)dt. 

respectively. 

3et 

Set  =  0. 

do  i  =  l,  M 
doj=l,p 

For  each  ,  compute 

)  .  [j,  +  vun-  ^(t)dt  +  »'«(<)  •  ^{i)di  +  J«] 

end  do 

Jl  =  Jl  +  idBr) 
end  do 


Thus,  the  amount  of  work  required  in  Step  6(a)  is  0{N  -n).  The  integrals  required  on  eacli 
subinterval  in  Step  6  (b)  can  be  computed  by  spectral  integration  (see,  for  example,  [ll])  using 
0{p-  logp  •  v)  work.  The  total  cost  is  therefore  of  the  order  0{M  -  p-  logp  ■  n)  or  0[i\  ■  logp  ■  n  )■ 

4.4.  Informal  Description  of  the  Algorithm  for  Nonlinear  DDEs 

The  nonlinear  algorithm  is  a  straightforward  application  of  the  Unear  algorithm  describeu  in 
the  previous  subsection.  The  solution  is  obtained  using  Newton’s  method  for  nonlinear  ODF.s; 
each  Newton  iterate  is  obtained  by  solving  the  Unearized  problem  (50)  via  the  algorithm  of  tlie 
preceding  subsection. 

As  with  .Mgorithm  A,  we  subdivide  the  interval  [a,c]  into  a  large  number  of  subinterval.s  M : 
for  simplicity  wc  assume  M  —  2"’,  with  m  a  positive  integer.  As  before  the  boundary  points 

^1-^2 . are  defined  by  (132),  and  the  intervals  B[,  [I  <  I  <  Jn),(l  <  i  <  2')  by 

(133). 

On  the  step  of  the  Newton  process.  Algorithm  A  is  applied  to  the  integral  equation 

P'‘h  =  9k,  (149) 

with  the  operator  :  (Z/^[a,c])”  — ♦  (Z,^[a,c])”  defined  by  the  formula 

P’'{Sk){x)  =  6k{x)  +  ilk(x)  ■  I  Go(xj)  ■  6k{t)dt,  (150) 
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with  Sk  the  solution  of  the  integral  equation  (50),  fit  given  by  (51),  and  gk  given  by  (52),  The 
integral  equations  (139),  (140)  now  assume  the  form 


(151) 

Qui4>i.l)  = 

(152) 

with  the  operator  Pj^i  :  (Z^[a,c])"  — ►  (jL^[a,c])”  defined  by  the  formula 

Pl‘/6k){x)  =  6k{x)  +  akiT)-  Go{x,i)-6k{t)dt, 

•'*!+(. -1)2'"-' 

(153) 

and  the  operator  Q'^i  :  (i^fa,  c])”’^"  — ►  (i,^[a, c])”’'”  defined  via  the  formula 

\(:r)  +  fi)t(a:)-  Go{x.t)-x(t)d1. 

•'*’l4(,-l)  2’"-' 

(151) 

Once  Algorithm  .4  has  computed  the  solution  6k  to  (149),  we  obtain  <Jk+i  (H)-  ^As-i 

via  (49). 

The  nonlinear  algorithm  requires  an  initial  approximation  $o-^o  to  the  solution  $  and  its 
derivative  4*'  of  equation  (5),  and  we  assume  that  both  are  supplied  by  the  calling  program, 
is  obtained  from  '"'S-  the  identity 

Mx)  =  ^o{a^)  +  Po(x)  •  ^o(i)- 

(1.55) 

The  procedure  is  terminated  when  the  stopping  criterion 

Ikl-lb  - 

(1.56) 

is  satisfied,  with  t  provided  by  the  calling  program.  Since  Newton’s  method  frequently  fails 
to  converge,  the  calling  program  also  permits  a  certain  maximum  number  of  iterations,  after 
which  the  algorithm  stops,  signalling  failure. 

The  following  is  a  more  detailed  description  of  the  .■  umerical  procedure. 

Algorithm  B 


Coniinent  [Define  the  computational  grid.] 

Create  M  =  2"'  subintervals  on  |o,c]  by  choosing  a  sequence  of  boundary  points  bi,bn . fc.i/.fc.ij  +  i 

with  bi  =  a  and  +  i  =  c.  Choose  the  number  p  of  Chebyshec  nodes  on  each  interval  =  [6,.  6,  +  )] 

for  j  =  1, . . . ,  A/.  Determine  the  locations  of  the  scaled  Chebyshev  nodes  t,’  ,  r,', . . , ,  on  each  interval 
and  use  the  initial  approximations  <I>o,4’o  fo  evaluate  the  initial  approximations  4>.<r  at  these  node 
obtaining  ,  cr,  ,n.  Choose  tolerance  c. 

Coniinent  [Use  Algorithm  A  to  compute  Newton  iterates  obtaining  the  solution  <l>  of  equation  (o).] 
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repeat 

(1)  Set  4>  =  $,  ff  =  cr. 

(2)  Evaluate  the  functions  Q,  g  at  each  of  the  scaled  Chebyshev  nodes  t}  ,Tf , . . .  on  each 
interval  ,  obtaining 

(3)  Apply  Algorithm  A  to  the  discretized  form  of  (50),  obtaining  6. 

(4)  Set  <T  =  O’  +  6. 

(5)  Compute  the  solution  4>  of  equation  (49)  from  the  values  of  <r,  by  using 
composite  Chebyshev  quadrature  (see  Step  6  of  Algorithm  A  and  Remark  4.5  above). 

until  ll^olb/lkolb  <  f 


V.  Numerical  Results 

FORTRAN  programs  have  been  written  implementing  both  algorithms  described  in  the 
preceding  section.  In  this  section,  we  discuss  several  details  of  our  implementation,  and  demon¬ 
strate  the  performance  of  tlie  scheme  with  numerical  examples. 

The  following  technical  details  of  our  implementation  appear  to  be  worth  mentioning. 

1.  The  algorithms  described  in  the  preceding  section  require  that  the  number  M  of  elementary 
subintervals  on  the  interval  [a,c]  be  a  power  of  2.  Clearly,  this  is  not  an  essential  limitation 
and  it  can  be  removed  by  simple  bookkeeping  changes.  In  the  version  of  the  algorithms  used 
for  numerical  experiments,  these  changes  were  made. 

2.  Algorithm  A  depends  for  its  stability  on  the  equations  (139),  (140)  having  unique  solutions 
for  all  subintervals  B[  (/  =  0,1,... ,M,  i  =  1,...,2^),  while  Algorithm  B  depends  on  (151). 
(152)  having  unique  solutions  for  all  subintervals  B\  and  for  all  Newton  iterates  k.  It  is  easy 
to  construct  examples  for  which  these  conditions  are  violated,  even  though  equation  (57)  or 
equation  (32)  has  a  unique  solution.  In  such  cases,  a  different  subdivision  of  the  interval  [u.c] 
can  be  attempted,  such  that  none  of  the  subintervals  B^  of  the  new  subdivision  coincides  witli 
an  interval  of  the  original  one.  This  procedure  can  be  viewed  as  a  form  of  pivoting,  and  it  i.s 
easy  to  show  that  it  is  always  possible  to  make  it  w’ork.  It  has  not  been  implemented  at  this 
point,  and  we  have  not  so  far  encountered  a  need  for  it. 

3.  We  have,  however,  implemented  a  crude  scheme  for  detecting  high  condition  numlrei.s  in 
the  algorithms.  These  can  occur  in  tw'o  places:  in  the  solution  of  the  linear  systems  on  eacli  of 
the  finest  level  subintervals  (Step  1  of  Algorithm  A),  and  while  computing  coefficients  A).  Ai 
defined  by  (90),  (91)  used  when  merging  solutions  on  two  consecutive  subintervals  (Step  3  of 
Algorithm  A).  In  both  cases,  the  condition  number  of  the  system  being  solved  is  estiinaled 
in  the  process  of  solution  (we  use  a  standard  LINPACK  routine),  and  the  largest  of  these  i.s 
returned  to  the  user.  When  an  extremely  large  condition  number  is  detected  by  the  LINP.ACK 
routine,  the  resulting  solution  of  the  original  ODE  should  be  viewed  as  suspect.  It  is  easy 
to  show  that  when  the  differential  operator  is  positive  definite,  this  cannot  happen.  A  more 
complete  treatment  of  this  subject  requires  further  study. 
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4.  In  the  upward  sweep  (Step  3)  of  Algorithm  A,  we  evaluate  the  matrices  for  all 

intervals  5,,;  and  use  these  matrices  to  evaluate  the  vectors  the  vectors  A*’^  and, 

finally,  the  solution  a  of  the  integral  equation  (57).  But  the  matrices  do  not  depend 

on  the  right-hand  side  /  of  equation  (57),  and  it  is  easy  to  see  that  their  evaluation  accounts 
for  more  than  90%  of  the  work.  Therefore,  whenever  the  equation  (57)  has  to  be  solved  with 
multiple  right-hand  sides,  we  can  precompute  the  matrices  and  store  them,  saving  90% 

of  the  cost  of  the  evaluation  of  subsequent  solutions. 

The  algorithms  of  this  paper  have  been  applied  to  a  variety  of  problems.  Seven  experiments 
are  described  below,  and  their  results  are  summarized  in  Tables  1-18. 

Tables  1-13  are  associated  with  examples  for  which  analytic  solutions  are  available.  In  eacli 
of  these  tables,  the  first  column  contains  the  total  number  N  of  nodes  in  the  discretization  of 
the  interval  [a,c].  The  second  column  contains  the  relative  error  of  the  numerical  solution 
as  compared  with  the  analytically  obtained  one  at  5000  equispaced  points  within  the  interval 
(a,  c],  where  Chebyshev  interpolation  has  been  used  to  evaluate  the  numerical  solution  at  eacli 
of  the  5000  points.  The  third  column  contains  the  maximum  absolute  error  obtained  at  any  of 
the  5000  points.  Columns  four  and  five  contain  the  same  information  for  the  derivative  of  the 
solution  (i.e.  its  relative  and  absolute  errors  respectively).  The  sixth  column  contains 
the  CPU  time  required  to  solve  the  problem,  excluding  the  time  used  to  evaluate  the  solution 
at  5000  equispaced  points,  where  in  all  cases  the  times  are  given  for  a  SUN  SPARCstation 
1  computer.  Tables  11-13,  associated  with  a  nonlinear  example,  have  in  addition  a  seventh 
column  which  contains  the  number  of  Newton  steps  taken  before  the  stopping  criterion  (156) 
has  been  satisfied,  with  e  =  10“'°. 

Tables  14-18  are  associated  with  examples  for  which  we  did  not  have  analytic  solutions.  In 
these  examples,  we  compare  each  numerical  solution  with  p  Chebyshev  nodes  and  n  subintei  vals 
against  the  solution  with  p  Chebyshev  nodes  and  2-n  subintervals.  In  each  of  these  tables,  the 
first  column  contains  the  total  number  N  of  nodes  in  the  discretization  of  the  interval  [a,  c].  Tlie 
second  column  contains  the  relative  error  of  the  numerical  solution  as  compared  with  the 
numerical  solution  with  twice  the  number  of  subintervals,  where  the  comparison  is  made  at  each 
of  5000  equispaced  points  in  the  interval  [a,c],  and  where  Chebyshev  interpolation  has  been 
used  to  evaluate  the  numerical  solution  at  each  of  the  5000  points.  The  third  column  contain> 
the  maximum  absolute  error  obtained  at  any  of  the  5000  points.  Columns  four  and  five  contain 
the  relative  and  absolute  Z*”  errors,  respectively,  for  the  derivative  of  the  solution.  The 
sixth  column  contains  the  SPARCstation  CPU  time  required  to  solve  the  problem,  excluding 
the  time  used  to  evaluate  the  solution  at  5000  equispace  points.  Tables  16-18,  associated  with 
a  nonlinear  example,  have  in  addition  a  seventh  column  which  contains  the  number  of  Newton 
steps  taken  before  the  stopping  criterion  (156)  has  been  satisfied,  with  c  =  10~'°. 

Remark  5.1  In  Examples  3-4  below,  we  solve  boundary  value  problems  of  order  2;  in  Exanijile 

5,  we  solve  a  problem  of  order  4;  and  in  Example  7  we  solve  a  system  consisting  of  four 
equations  of  order  2  and  two  equations  of  order  4.  In  all  these  cases,  the  problems  were 
reduced  to  canonical  first  order  systems  (see,  for  example,  [4]),  with  the  latter  solved  by  means 
of  algorithms  A  or  B  of  the  preceding  section,  as  appropriate. 
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n 

£■"($) 

£«>($) 

E~($') 

t  (sec.) 

16 

0.962  xlO" 

0.685  x10* 

0.261  xlO* 

0.673  XlO'*  " 

0.150  xlO^ 

32 

0.244  xlO^ 

0.216  xlO^ 

0.108  xlO^ 

0.169  xlO® 

0.300  xl0° 

64 

0.700  XlO-^ 

0.272  XlO^ 

0.200  XlO^ 

0.261  XlO'* 

0.560  xl0° 

128 

0.255  XlO-i 

0.125  xlO* 

0.754  Xl0° 

0.125  XlO” 

0.108  xlO^ 

256 

0.667  xlO-2 

0.342  Xl0° 

0.198  Xl0° 

0.342  xlO^ 

0.214  xlO^ 

512 

0.805  XlO"^ 

0.536  xlO'^ 

0.238  XlO"^ 

0.536  XlO^ 

0.428  xl0‘ 

1024 

0.330  XIO-^ 

0.291  xlO-2 

0.977  xlO-3 

0.291  XlO' 

0.846  xlO’ 

2048 

0.482  X10-® 

0.529  xlO-"* 

0.143  XlO-^ 

0.529  xlO"’ 

0.168  xlO^ 

4096 

0.316  xlO-* 

0.476  X10-® 

0.935  xlO-^ 

0.476  xlO-^ 

0.337  xl02 

8192 

0.112  xl0-^° 

0.170  xlO-® 

0.331  X10-® 

0.170  xlO-5 

0.670  xlO^ 

16384 

0.115  X10-” 

0.113  xl0-’° 

0.264  X  10-^1 

0.140  xlO-^ 

0.137  xlO^ 

Table  1:  Numerical  results  for  Example  1,  p  =  S. 

Example  1  This  example  is  taken  from  [2],  where  it  is  introduced  as  a  stiff  problem.  The 


equation  to  be  solved  is  given  by  the  formulae 

998- 1998- (Ji2(j)  =  2-x,  (3  57) 

<;l)2(i)  +  999  •  +  1999  •  <;i'2(x)  =  x,  (loS) 

subject  to  the  boundary  conditions 

-^i(O)  =  1,  (159) 

4>2(1)  =  -6-e-' +5-e“^°°®  +  .004-(.999  +  .00]  (IGO) 


We  use  the  results  of  Theorem  2.3  to  reduce  the  first  order  system  (157)  -  (160)  to  one  subject 
to  homogeneous  boundary  conditions 


0i(O)  =  0,  (IGl) 

d>2(l)  =  0.  (1G2) 

We  apply  Algorithm  A  to  this  system  using  equispaced  subintervals,  with  the  number  of  Chein  - 
shev  nodes  p  =  8,16,24.  For  this  experiment,  the  background  Green’s  function  is  chosen  to 
correspond  to  the  equation 

$'(x)  =  0,  (1G3) 

subject  to  boundary  conditions  (161)-(162).  The  results  of  this  experiment  are  presented  in 
Tables  1-3. 

Example  2  We  solve  the  problem  (157)  -  (160)  defined  in  Example  1,  but  using  an  alternate 
division  of  the  subintervals.  Since  the  solution  of  this  problem  has  a  fairly  sharp  boundary 
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n 

W(J) 

£oo($) 

E^) 

16 

0.567  x10" 

0.436  xlQi 

0.916  Xl0° 

0.397  xlO'i 

0.310  xio" 

32 

0.251  xlO^ 

0.192  XlO^ 

0.673  xlOi 

0.175  xlO® 

0.540  XlO" 

64 

0.340  xlO-* 

0.139  xlQi 

0.859  xlO" 

0.129  xlO'’ 

0.950  XIO" 

128 

0.744  xlO-2 

0.344  Xl0° 

0.220  XlO" 

0.344  xl03 

0.179  xlQi 

256 

0.798  xlO-3 

0.389  xlO-i 

0.236  xlO-i 

0.389  xl02 

0.345  XlO’ 

512 

0.164  xlO-'* 

0.930  xlO-2 

0.486  xlO-3 

0.930  xl0° 

0.686  xlOi 

1024 

0.364  XIO"^ 

0.243  xlO-3 

0.108  xlO-3 

0.243  xlO-2 

0.137  xl02 

2048 

0.992  X10-" 

0.817  xlO-" 

0.294  xlO-® 

0.817  xlO-® 

0.274  Xl02 

4096 

0.942  X 10-13 

0.776  X 10-12 

0.146  xlO-12 

0.753  xlO-" 

0.532  Xl02 

8192 

0.214  xlO-13 

0.164  XlO-ii 

0.328  xlO-12 

0.164  xlO-® 

0.107  Xl03 

Table  2:  Numerical  results  for  Example  1,  p  =  16. 


71 

■Earn 

t  (sec.) 

24 

0.511  xlO" 

0.373  xlOi 

0.752  XIO" 

0.358  xlO"! 

0.690  xlO" 

48 

0.251  xlO" 

0.244  XlOi 

0.112  xlQi 

0.200  xlOi 

0.116  xlQi 

96 

0.591  xlO-2 

0.240  xlO° 

0.175  xlO" 

0.240  xl03 

0.209  xlOi 

192 

0.496  xlO-3 

0.214  xlO-i 

0.147  xlO-i 

0.214  Xl02 

0.387  xlOi 

384 

0.546  xlO-3 

0.258  xlO-3 

0.162  xlO-3 

0.258  xlO" 

0.765  xlO’ 

768 

0.274  xlO-® 

0.129  xlO-" 

0.811  xlO-^ 

0.129  xlO-3 

0.147  Xl02 

1536 

0.105  xlO-12 

0.335  xlO-ii 

0.142  xlO-ii 

0.256  xlO-® 

0.295  Xl02 

3072 

0.663  X 10-13 

0.624  X 10-12 

0.125  xlO-12 

0.621  XlO-" 

0.578  Xl02 

Table  3:  Numerical  results  for  Example  1,  p  =  24. 
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n 

W{¥) 

£“($) 

W(W) 

t  (sec.) 

16 

0.567  xl0“ 

0.436  xl0‘ 

0.916  xlO" 

0.397  xlO" 

0.310  xlO^ 

32 

0.251  xlO^ 

0.192  Xl02 

0.673  xlOi 

0.175  xlO^ 

0.540  xlO" 

64 

0.790  XlO-2 

0.360  Xl0° 

0.222  xl0“ 

0.345  xlO^ 

0.950  Xl0° 

128 

0.992  X10-” 

0.818  xlO-3 

0.294  xlO"® 

0.816  X10-® 

0.177  xlOi 

256 

0.244  X  10-^2 

0.261  xlO-ii 

0.525  X 10-12 

0.228  xlO-® 

0.352  xlOi 

512 

0.243  X 10-12 

0.261  xlO-" 

0.525  xlO-12 

0.229  xlO-® 

0.681  xlOi 

Table  4:  Numerical  results  for  Example  2,  p  =  16. 


n 

W(¥j 

£“($) 

F(¥) 

t  (sec.) 

24 

0.511  Xl0“ 

0.373  xlOi 

0.752  Xl0° 

0.358  xlO"! 

0.710  xlOii 

48 

0.251  xl0° 

0.244  xlOi 

0.112  xlQi 

0.200  xlO"! 

0.117  xlOi 

96 

0.496  xlO-3 

0.214  xlO-i 

0.147  xlO-i 

0.214  xl02 

0.209  xlO’ 

192 

0.293  X 10-12 

0.227  xlO-ii 

0.507  X 10-12 

0.229  XIO"® 

384 

0.294  X  10-12 

0.227  xlO-ii 

0.509  X 10-12 

0.230  xlO-® 

Table  5:  Numerical  results  for  Example  2,  p  =  24. 

layer  near  the  left  end  of  the  interval  [0,1],  we  construct  the  intervals  =  [6,,  6,+]]  via  ihc 
formula 


6i  =  0, 

/  Af+l-ti 

6,  =  (-]  for  t  =  2,...,M+ 1  ,  (164) 

so  that  they  become  progressively  smaller  near  the  left  end  of  the  interval  [0.  ij.  As  in  Example 
1,  we  reduce  the  problem  (157)-  (160)  to  a  first  order  system  subject  to  homogeneous  l)ound- 
ary  conditions  (161)-(162).  Algorithm  A  has  been  applied  to  this  problem  using  the  Gieen's 
function  corresponding  to  the  equation  (163)  subject  to  boundary  conditions  (161)-(162).  and 
with  the  number  of  Chebyshev  nodes  p  =  16  and  24.  The  results  of  this  c.xpcrimeiit  appe.ir  in 
Tables  4-5.  and  are  most  satisfactory. 

Example  3  This  example  is  taken  from  [15],  where  it  is  described  as  a  reasonably  difficult  oiu 
due  to  the  presence  of  rapidly  growing  solutions  of  the  corresponding  homogeneous  equation. 
The  equation  to  be  solved  is 

(p”  -b  400  ■  4>  =  -400  •  cos^(7r  •  x)  —  2  ■  ■  cos(2  •  tt  ■  x),  ( 1(>5) 

subject  to  the  boundary  conditions 

0(0)  =  <^(1)  =  0.  (166) 
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Table  6:  Numerical  results  for  Example  3,  p  =  S. 


Table  7:  Numerical  results  for  Example  3,  p  =  16. 

We  reduce  the  problem  (165),  (166)  to  a  first  order  system,  and  apply  Algorithm  A  to  this 
system  using  equispaced  subintervals,  with  the  number  of  Chebyshev  nodes  p  =  8, 16  and  24. 
For  this  experiment,  the  background  Green’s  function  is  chosen  to  correspond  to  the  equation 

^’'(x)+(_^l  ■^^•$(x)  =  0.  (167) 

subject  to  boundary  conditions  (166).  The  results  of  this  experiment  are  presented  in  Tal)les 
6-8. 

Example  4  We  solve  the  problem  (165),  (166)  defined  in  Example  3,  but  we  use  the  back¬ 
ground  Green’s  function  corresponding  to  the  equation 

4>'(i)  =  0,  (168) 

subject  to  boundary  conditions  (166).  We  reduce  the  problem  (165),  (166)  to  a  first  order 
system,  and  then  use  the  results  of  Theorem  2.9  to  express  by  the  formula 

$(i)=  «'(i).r(i). 
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n 

I^(¥) 

W{¥) 

E~(4>') 

t  (sec.) 

24 

48 

96 

0.164  XIO"^ 
0.338  Xl0-i‘‘ 
0.264  XIO-^'* 

0.336  xlO-^ 
0.568  xlO-13 
0.853  X  10-^3 

0.156  xlO-« 
0.345  XlO-^'' 
0.273x10-^“ 

0.676  XlO-^ 
0.188  xlO-" 
0.261  xlO-" 

0.760  xlO^ 
0.136  xlO’ 
0.257  xlO’ 

Table  8:  Numerical  results  for  Example  3,  p  =  24. 


with  '5’  ;  [0, 1]  — ♦  L(R^’‘^)  given  by  the  formula 

«'(x)  = 


1  —  X  X 
-1  1 


and  r  :  [a,c]  — >  R”  the  solution  to  the  equation 

r'(i)  +  .  ('lr'(i)  +  p{x)  ■  <f>{x))  ■  r(x)  =  ■  fix), 

subject  to  boundary  conditions 


1  0 
0  0 


^'(a)  •  r(a)  + 


0  0 
0  1 


’j'(c)-r(c)  =  o, 


with  the  matrix  valued  function  p  defined  by  the  formula 


p{x)  = 


0  -1 


400  0  J 

and  the  vector  valued  function  f  defined  via  the  formula 

/ 

0 


i  -400  •  cos^(7r  •  i)  -  2  •  •  cos(2  ■  tt  •  x) 


(169) 


(170) 


[171' 


(1721 


[1731 


We  apply  Algorithm  A  to  this  problem  using  equispaced  subinlervals,  with  the  number  of 
Chebyshev  nodes  p  =  16  and  24.  The  results  of  this  experiment  are  presented  in  Tables  9-10. 


Example  5  We  consider  a  system  of  Jacobian  elliptic  functions  sn,cn,(In  :  [0,10  •  K]  —  R 
(see,  for  example,  [l])  which  are  solutions  to  the  equations 


sn'(i)  =  cn(x)  •  dn(x),  (174) 

cn'(i)  =  —sn{x)  ■  dn(x),  (17')) 

dn'{x)  =  — m  •  sn(i)  •  cri(x),  (176) 

with  m  =  1  in  our  experiments,  subject  to  the  boundary  conditions 

5n(0)  =  0,  (177) 

cn{0)  =  1,  (178) 

cfTi(40-/v')  =  1,  (179) 
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n 

E^) 

t  (sec.) 

16 

0.946  XlO-^ 

0.193  XlO-2 

0.887  xlO-'* 

0.412  XIO"^ 

0.320  xlO'i 

32 

0.281  xlO-^ 

0.567  xlO-® 

0.307  xlO-^ 

0.145  xlO-'i 

0.530  xl0° 

64 

0.517  xlO-” 

0.833  xl0-*° 

0.527  xlO-” 

0.198  X10-® 

0.960  X10° 

128 

0.106  xl0-^‘‘ 

0.355  X 10-13 

0.139  xlO-i*! 

0.796  X 10-13 

0.176  xlQi 

256 

0.977  xlO-^® 

0.391  xlO-13 

0.140  xlO-i’* 

0.853  X 10-13 

0.347  xlQi 

Table  9;  Numerical  results  for  Example  4,  p  =  16. 


n 

f;3($) 

E°°{^) 

W(¥) 

£;«($/) 

t  (sec.) 

24 

48 

96 

0.113  xlO-® 
0.470  Xl0-i‘‘ 
0.180  xlO-i-* 

0.232  xlO-*^ 
0.711  xlO-13 
0.391  X 10-13 

0.106  xlO-® 
0.371  xlO-i'i 
0.187x10-1“ 

0.484  XlO-*^ 
0.176  xlO-ii 
0.909  X 10-13 

0.700  xlQii 
0.116  xlOi 
0.207  xlOi 

Table  10:  Numerical  results  for  Example  4,  p  =  24. 


with  K  given  by  the  expression 


de 


\/l  —  m  •  sin^  9 


;]so) 


We  use  for  an  initial  guess  the  solution  to  (176),  (179)  for  m  =  0,  which  is  defined  by  tlio 
formulae 


sn{x)  =  sin  . 

cn(i)  = 


dn(x)  =  1. 


(151) 

(152) 

(153) 


We  .reduce  the  problem  (165),  (166)  to  a  first  order  system,  and  then  use  the  result.^  of 
Theorem  2.4  to  reduce  this  system  to  one  subject  to  the  homogeneous  boundary  condition' 


5n(0)  =  0,  (181) 

cn(0)  =  0, 

dn{40K)  =  0.  (186) 

We  then  apply  Algorithm  B  to  this  system  using  equispaced  subintervals,  with  the  number  of 
Chebyshev  nodes  p  =  8, 16  and  32.  For  this  experiment,  the  background  Green’s  function  is 
chosen  to  correspond  to  the  equation 

$'(z)  =  0, 
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n 


n 

E'\^) 

f;°°($) 

E°^(<f>') 

t  (sec.) 

Ste])s 

128 

0.211  xlO-^ 

0.551  xlO'^ 

0.333  xlO-^ 

0.620  X10-’ 

0.167  xlO'^ 

9 

256 

0.158  xlO-2 

0.416  xlO-2 

0.228  xlO-2 

0.430  xlO-2 

0.258  X102 

7 

512 

0.849  X10-® 

0.226  xlO-^ 

0.123  XlO-** 

0.240  XIO-" 

0.439  XlO^ 

6 

1024 

0.106  xlO-’’ 

0.330  xlO-^ 

0.180x10-'^ 

0.460  XIO"^ 

0.883  xlO^ 

6 

2048 

0.313  xl0-^° 

0.984  xl0-^° 

0.599  Xl0-*° 

0.150  xlO"® 

0.175  XIO^ 

6 

4096 

0.147  xlO-^2 

0.469  Xl0-'2 

0.264  X  10-^2 

0.674  xl0-'2 

0.348  xlO^ 

6 

Table  11:  Numerical  results  for  Example  5,  p  =  8. 


73 

E^{^) 

WW) 

t  (sec.) 

Steps 

64 

0.162  X10“ 

0.505  Xl0“ 

0.244  X10° 

0.546  X10° 

liiiirMmai 

10 

128 

0.422  xlO"^ 

0.108  xl0° 

0.612  X10-* 

0.113  Xl0° 

0.246  xlO^ 

i 

256 

0.181  XIO"^ 

0.476  xlO-3 

0.268  xl0-° 

0.491  xl0-° 

0.409  xlO^ 

6 

512 

0.441  xl0-‘ 

0.120  xlO-® 

0.672  xlO-"^ 

0.159  xl0-° 

0.851  xlO^ 

6 

1024 

0.425  XlO"'^ 

0.125  xlO'” 

0.115  xlO-^i 

0.293  X10-” 

0.164  xl0° 

6 

2048 

0.115  xlO-^2 

0.317  xl0-’2 

0.164  XlO-12 

0.317  xlO-^2 

6 

4096 

0.569  xlO-*^ 

0.152  xl0-'2 

0.814  Xl0-^° 

0.154  xl0-'2 

6 

Table  12:  Numerical  results  for  Example  5,  p  =  16. 


subject  to  boundary  conditions  (184)  -  (186).  The  results  of  this  experiment  are  presented  in 
Tables  11-13. 

Example  6  This  example  is  taken  from  [14].  Its  purpose  is  to  demonstrate  the  performance 
of  the  method  when  the  equation  to  be  solved  contains  fourth  order  derivatives.  The  deflection 
of  a  beam  under  a  uniform  load  q,  with  the  beam  built  in  at  the  left  end  (x  =  D)  and  sinii)!v 
supported  at  the  right  end,  is  given  by  the  formula 


s"'V)+A-jn)=  ’ 


E  l  E  l' 

subject  to  the  boundary  conditions 

y{0)  =  y'{0)  =  y{L)  =  y"(L)  =  0, 


(  l.SS) 


with  k  the  force  per  unit  deflection  per  unit  length  of  beam,  and  E  ■  I  the  flexural  rigidity  of 
the  beam.  The  values  of  the  constants  used  are 

=  1.2  X  10^  in., 

=  2.604  X  10^  psi. 


L 

k 
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(18t)) 

(190) 


n 

£~($') 

t  (sec.) 

Steps 

256 

6 

512 

B  i  &  M 

6 

1024 

QBQRBS 

B  ^3  3 

6 

2048 

6 

Table  13;  Numerical  results  for  Example  5,  p  =  32. 


4.34  X  10^  Ibs/in., 

(191) 

3.0  X  10^  psi. 

(192) 

3.0  X  10^  in.^ 

(193) 

(see  [14],  p.  174).  The  norm  of  y,  the  solution  to  (188),  is  approximately  10®  times  larger 
than  tlie  norm  of  y"".  Combined  with  the  high  number  of  derivatives  in  (187),  this  tends  to 
present  difficulties  for  finite  difference  methods.  We  reduce  the  problem  (187),  (188)  to  a  first 
order  system,  and  then  use  the  results  of  Theorem  2.9  to  express  $  by  the  formula 


$(l 

H 

II 

r(x 

), 

with  ^(i) 

:[0, 

i]-< 

L(R‘‘ 

*‘‘‘)  given  by  the  formula 

/  I'-r 
L 

0 

X 

L 

0 

\ 

^(x) 

0 

1 

0 

0 

— 

0 

0 

0 

1 

1  -1 

0 

1 

0 

/ 

and  r  :  [o, 

R" 

the  solution 

to  the  equation 

r'( 

i)  +  'J' 

-i(x) 

x) +  p(x) 

■  n 

a^))- 

r( 

subject  to 

boundary  conditions 

/  1 

0  0 

0  \ 

1 

/  0 

0 

0 

0 

0 

1  0 

0 

^(a) 

•r( 

«)  + 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

1 

0 

1  0 

0  0 

0 

V  1 

0 

0 

0 

(194; 


195 


^(c)  ■  r(c)  =  0. 


(196) 


with  the  matrix  valued  function  p  defined  by  the  formula 

p(i)  = 


f  ^ 

-1 

0 

0 

\ 

0 

0 

-1 

0 

0 

0 

0 

-1 

k 

\  FI 

0 

0 

0 

y 

(197) 
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n 

eW) 

t  (sec.) 

4 

0.115  XlO*^ 

0.672  xl0-‘ 

0.277  xl0-‘ 

0.433  xlO-^ 

0.110  xl0° 

8 

0.755  xlO-2 

0.506  xlO-2 

0.215  xlO-2 

0.426  XlO-"* 

0.220  Xl0° 

16 

0.474  xlO-3 

0.369  xlO-^ 

0.141  XlO-3 

0.317  xlO-® 

0.400  Xl0° 

32 

0.269  xlO-'* 

0.249  xlO-'* 

0.891  xlO-5 

0.212  X10-® 

0.780  xl0° 

64 

0.185  xlO-® 

0.162  xlO-^ 

0.554  X10-® 

0.137x10“'^ 

0.153  xlO' 

128 

0.117  xlO-® 

0.103  xlO"® 

0.338  XIO""^ 

0.866x10-® 

0.303  xlO' 

256 

0.891  xlO-* 

0.652  xlO-® 

0.223  X10-® 

0.539x10-^° 

0.603  xlO’ 

512 

0.306  XIO"® 

0.170  xlO"® 

0.338  X10-® 

0.416x10"'° 

0.120  xlO^ 

Table  14;  Numerical  results  for  Example  6,  p  =  4. 


n 

t  (sec.) 

8 

0.130  xlO"^ 

0.816  xlO"^ 

0.688  xl0-° 

0.139  XlO-' 

0.280  xl0° 

16 

0.792  xlO"® 

0.510  xlO"® 

0.101  xio-'^ 

0.149  XlO"® 

0.530  xl0° 

32 

0.465  xlO"® 

0.235  xlO"® 

0.468  XlO"® 

0.565  xl0-'° 

0.105  xlO' 

64 

0.916  xlO"® 

0.463  xlO"® 

0.927  X10-® 

—  1 

0.111x10-® 

0.201  XlO’ 

Table  15:  Numerical  results  for  Example  6,  p  =  S. 


and  the  vector  valued  function  /  defined  via  the  formula 


/  0  \ 
0 
0 

V  ^  I 


(lbs) 


We  apply  Algorithm  A  to  this  problem  using  equispace  subintervals,  with  the  number  of 
Chebysliev  nodes  p  =  4  and  8.  For  this  experiment  the  background  Green’s  function  is  chosen 
to  correspond  to  the  equation 

$'(x)  =  0, 

sub  ject  to  boundary  conditions  (188).  The  results  of  this  experiment  are  presented  in  'lables 
14-15. 


i 


Example  7  This  example  is  taken  from  [13],  and  was  developed  in  cooperation  with  the  author 
of  [13].  The  problem  is  a  system  of  six  nonlinear  equations  with  inhomogeneous  boundaiy 
conditions,  and  is  described  in  detail  in  Appendix  A.  Since  2  of  the  ODEs  are  fourth  order, 
and  4  are  second  order,  the  problem  reduces  to  a  system  of  16  first  order  nonlinear  equations. 
The  2  fourth  order  ODEs  have  in  their  fourth  derivatives  a  boundary  layer  of  order  lO'’  at  the 
left  end  of  the  interval. 
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n 

F(F) 

£0O(^/) 

t  {sec.) 

Steps 

80 

0.697  xlO"^ 

0.516  xlO^ 

0.118x10° 

0.415  xlO'* 

0.249  XIO^ 

6 

160 

0.854  xlO-2 

0.773  xlQJ 

0.167  xlO-i 

0.790  xlO^ 

0.542  xlO^ 

6 

320 

0.795  xlO-^ 

0.952  Xl0° 

0.177  xlO-2 

0.111  xlO^ 

0.105  XIO'* 

6 

640 

0.589  xl0-‘' 

0.917  xlO"^ 

0.143  xlO-3 

0.114  xl02 

0.205  xlO"* 

6 

Table  16:  Numerical  results  for  Example  7,  p  —  4. 


n 

WiJ) 

£oo($) 

t  (sec.) 

Steps 

80 

0.814  XlO-'-^ 

0.686  xlO* 

0.151x10-’ 

0.667  xlO^ 

0.633  xlO'^ 

6 

160 

0.409  XlO-2 

0.380  xlO’ 

0.871  xlO-3 

0.454  xlO^ 

0.107  xlO'’ 

6 

320 

0.884  xlO-5 

0.977  xlO-2 

0.238  xlO-'* 

0.140  xio’ 

0.206  xlO' 

6 

640 

0.930  X10-" 

0.130  XIO-^ 

0.312  xI0-° 

0.213  XIO"’ 

0.396  XlO-* 

6 

Table  17:  Numerical  results  for  Example  7,  p  =  8. 

We  use  for  an  initial  guess  the  solution  to  the  problem  for  n  =  40, p  =  4.  We  apply 
Algorithm  B  to  this  system  using  equispaced  subinter%’als,  with  the  number  of  Cliebvsliev 
nodes  p  =  4,8,  and  IG.  The  results  of  this  experiment  are  presented  in  Tables  IG-IS. 

The  following  observations  can  be  made  from  Tables  1  -  18,  and  are  corroborated  by  our 
more  extensive  experiments. 

1.  The  practical  convergence  rate  of  the  method  is  consistent  with  the  theoretical  one.  f  or 
larger  p,  the  exact  numerical  verification  of  the  order  of  convergence  tends  to  be  difficult,  since 
the  precision  of  calculations  is  exhausted  before  the  behavior  of  the  scheme  becomes  asymptotic. 
However,  this  is  often  encountered  when  dealing  with  rapidly  convergent  algorithms. 

2.  For  small-scale  problems  and  large  p,  the  algorithm  produces  essentially  exact  results  with 
a  small  number  of  nodes.  For  large-scale  problems,  double  precision  accuracy  is  achieved  at 
approximately  20  nodes  per  wavelength  with  p  =  20,  at  12  nodes  per  wavelength  with  p  =  24, 


71 

£00(^) 

E^{7i>') 

t  (sec.) 

80 

0.673x10-3 

0.155  xlO"* 

6 

160 

0.291  XlO* 

6 

0.565  XlO"* 

6 

B 

0.112  xlO^ 

6 

Table  18:  Numerical  results  for  Example  7,  p  =  16. 
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and  at  10  nodes  per  wavelength  with  p  =  32.  The  optimal  timings  are  achieved  at  p  1)01  ween 
24  and  32  (provided  that  about  10  -  12  digits  of  accuracy  are  desired). 

3.  The  condition  number  of  a  Nystrdm  discretization  of  a  second  kind  integral  equation  is 
asymptotically  bounded,  and  our  results  reflect  this  fact.  The  relatively  poor  accuracy  (S  - 
11  digits)  obtained  in  Examples  6  and  7  is  due  to  the  ill-conditioning  of  the  original  ODE.  as 
opposed  to  that  of  the  numerical  scheme  used. 

4.  The  algorithm  is  completely  indifferent  to  the  stiffness  near  the  left  end  of  the  interval  [0. 1] 
of  equations  (157),  (158)  in  Examples  1-2. 

5.  It  is  easy  to  use  the  algorithm  in  an  adaptive  manner,  as  demonstrated  in  Example  2. 
However,  a  fully  adaptive  version  of  the  scheme  has  not  been  implemented.  The  intervals  Bl'‘ 
in  ExampE  2  were  provided  by  the  calling  program  (as  opposed  to  having  been  constructed  b}’ 
the  algoritnm  itself). 

6.  The  numerical  advantages  of  one  background  Green’s  function  over  another  tend  to  be  minor, 

as  indicated  in  Examples  3-4.  However,  using  the  Green’s  function  given  in  Lemma  2.2  doe.- 
result  in  a  slightly  faster  algorithm.  This  is  because  this  Green’s  function  is  constant  in  each 
of  the  intervals  (j  <  i),(x  >  i),  which  provides  in  Step  2  of  Algorithm  A  faster  evaluation  of 
the  matrices  given  by  equations  (142)-(143)  and  vectors  (!)[’''*,  6^'”  given  by  ecpiations 

(144)  (145).  and  provides  in  Step  6  a  faster  evaluation  of  the  solution  ^  of  equaiioji  (27). 

7.  The  algorithm  can  solve  systems  of  high  order  equations  with  no  numerical  difficult^-,  as 
demonstrated  by  Examples  6-7. 


VI.  Conclusions 

An  algorithm  has  been  presented  for  the  solution  of  two-point  boundary  value  problom.s  of 
ordinary  differential  equations.  The  algorithm  is  based  on  reducing  the  differential  equation  to 
a  second  kind  integral  equation,  with  the  subsequent  solution  of  the  latter  via  a  Nystrom  tyjro 
scheme.  It  has  CPU  time  requirements  proportional  to  N-p^-n^,  with  A’  the  number  of  node.s  in 
the  discretization  of  the  interval  of  definition  of  the  equation,  p  the  desired  order  of  convergenn' 
of  the  scheme,  and  n  the  number  of  equations  in  the  first  order  system.  The  method  doe.s  not 
iin  ohc  the  solution  of  linear  systems  with  large  condition  numbers,  permits  the  use  of  schemes 
with  extremely  high  orders  of  convergence,  and  is  quite  insensitive  to  boundary  layers  or  to 
end-point  singularities  in  the  coefficients  of  the  differential  equation. 

The  algorithm  has  been  combined  with  Newton’s  method,  resulting  in  a  sclieme  for  the 
solution  of  boundary  value  problems  for  nonlinear  ODEs.  In  this  case,  each  Newton  iterate  is 
exjji'ossed  as  the  solution  of  a  linear  second  kind  integral  equation;  the  analytical  and  numerira! 
advantages  of  integral  equations  are  thus  obtained  for  nonlinear  boundary  value  problem.'-. 
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Appendix  A. 

Numerical  Solution  of  Square-Cell  Convection  with  Strongly  Variable  Viscosity 


This  problem  is  introduced  in  [13]  as  a  model  for  square-cell  convection  with  strongly 
variable  viscosity.  The  model  is  used  to  examine  the  influence  of  temperature-dependent  and 
pressure-dependent  viscosity  on  convective  heat  transfer  and  surface  motion.  The  problem  is 
given  by  the  system  of  six  coupled,  nonlinear  equations 

=  (4•/^•0l(I)-|-4-(2•a^<(^)-a^,A,(x)-l-7.(<^"(a■).r'(x)  (190) 

-I-  2  ■  4>Tix) .  T'{x)  -2-al.  4>\(x)Jr\x)  -  2  ■  al  ■  <^';(x)  •  T(x) 

+  al  .  T"(x)  •  4>,{x)  -f-  at  ■  d>,(x)  •  T(i)))  -|-  2  •  7  •  i4>2{x)  •  e'^x) 

+  2  ■  d>2'(x)  •  -t-  4>2"{x)  ■  6i{x)  ~  2  ■  al  ■  4>2(x)  ■  6[(x)  -  3  ■  al  ■  o'^ix)  ■  0]{x) 

+  2-al-  4)2{x)  ■  e"{x)  +  2 -at-  <f>2{x)  ■  0i(x)  -2-al- 

-2-al-  ■  02(x))  -f  7  •  ■  e'^{x)  -H  2  •  •  ^2(^)  +  v'"[x)  -  diix) 

-3-al-  t'^'(x)  •  e2{x)  -  A-al-e'2{x)  -  tA(x)))/(4  •  (1  -  1 -T{x))). 

=  (2  •  +  2  •  (4  •  •  <j!)2(x)  -  4  •  •  <ji2(a-) -f  7  •  (^^'(x)  ■  r”(x)  (200) 

+  2  •  o'i'[x)  ■  r'(x)  -4-al-  <>'2(x)  •  T'(x)  -  4  ■  u|.  •  d>''(x)  •  T(x) 

+  2  •  al  -  T'\x)  -  d)2(x)  +  4  •  at-<?2{x)-T(x)  -b  d>'i'(x)  •  e'{{x)  -b  2  •  o'/'(.r )  •  0\{x) 

+  ■  ^1(2-')  -4-al-  •  ^'1(2;)  -  3  •  Ufc  •  d>'/(i)  •  ^i(x) 

+  al-<t>i(x).  <(x)  -b  2  •  Mx)  ■  0iix)))+  7  •  (- v'(x)  • 

-  2  ■  v"(x)  •  e'lix)  -  Xl^"'{x)  -  ei{x)  4-  4  -  al  ■  xl^'{x)  -  0i{x)  +  d  ■  al  -  0[[x)  -  xix))) 
/(2.(l-7-T(x))), 

r"(x)  =  al-{<i)[(x}-0i{x}-j-<^>i(x}-e[(x)4-^-(t>2ix)-02(x)4-^-62i-r)-02{x)).  (201) 

„2  _ 

^"(•2^)  =  {0i{x)  +  r'(x)-4>i{x))  +  2-e2{x)-(:>\{x)  +  2-0'2(x)-o,{x}  (202,. 

+  2  •  <?'2(X)  •  0i(x}  +  4-0[{x)-  <^2{x)  -  V>(x)  •  02(2-)). 

02(x}  =  al  ■{2-02(x}  +  2-T\x)-<^2(i)  +  2-0[{x)-(t>^(x)  +  Xix)-0i{x)).  (20.1) 

=  h  ■(4>i{x)-e'2{x)  +  <!>'l'{x}-02(x)-3-al-<l>\ix)-02{x)  +  al-<Pi[x)-0'2{x)  (201) 

-  <l>'^{x)  -  0[ix)  -  <t>'^'{x)  ■  e^{x)  4- 4-al-  <^2(2^)  ■0i{x)-2-  al  ■  <t>2[x)  -  ^[(x)) 

-b  10  •  (5  •  iix)  4-  7  •  (^'(x)  -  T'(x)  -  5  -  al  -  tix)  ■  fix))} 

+  J  •  (-4  •  ti'(2:)  •  eijx)  +  7-al-  Hx)  -  92{x)))/il0  -  10  ■  7  •  f(x)  -b  2  •  7  •  <  )). 
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subject  to  boundary  conditions 

Ma)  =  4>'l{a)  =  Ma)  =  </>2(a)  =  r(a)  -  ^  =  e,{a)  =  02(g)  =  ^'(a) 

Me)  =  4>'l{c)  =  d>2(c)  =  Me)  =  ^(c)  _  i  =  0i(c)  =  02(c)  = 

and  with  x  G  [—1/2, 1/2].  We  use  as  an  initial  guess  the  approximations 


4>\{x)  =  (  •  cos(7r  •  ar), 


(207) 


M^)  =  •  C2  •  sin(2  •  JT  •  x), 


(208) 


T(x)  =  •  C3  •  sin(2  •  TT  •  x)  -  I, 


0i(i)  =  f  •  Cl  •  cos(7r  •  x), 


02(x)  =  •  C4  •  sin(2  •  TT  •  x), 


3  ,  cos(7r  •  x)  +  jDd  •  cosh(\/5  •  a*:  •  x) 

=  c  •  7  •  vCs  - - :: - : - - 


TT^  +  5  •  a\ 


cos(3  •  -K  ’  x)  -  Z  ■  Dd  ■  cosh(  \/lak^) , ,,,  2, 


9'k‘^  +  5a? 


•)/(5g(.). 


The  constants  in  this  experiment  are  given  by  the  formulae 

R  =  2000, 

Gfc  =  2.2, 

Tr  =  10, 

.  -  2-(r.-l) 

Cr  +  1 

Too  = 


(213) 

(214) 

(215) 


=  "2  2' 

7r2  +  a2 


C2  = 


—  Ti-a^-C\  ■  Too 
2  •  ■  Too  -  8  •  Too  •  al ' 


at  ■  Cl 

C3  -  -7 - , 

4  ■  7r 

al-(2-C2  +  Tr -Cl) 
e^  —  n  2  ’ 

2  •  Too  •  af. 
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(222) 


1 


t 


i 


2 

Cs  =  .(y  •(7r’-5-a|)-Ci  •C2-(2-ir^  +  5-a^)), 

C6  =  TT  •  y  •  (y  •  (3  •  TT^  +  Cfc)  -  3  •  Cl  •  C2  •  Too  •  a|), 

-Dd  =  7r/(\/5  •  ajc  •  sinh(\/5  •  afc/2)), 

R  -  ^00 

TT  •  (C3  +  C4/4)  •  Too’ 

In  order  to  reduce  the  problem  to  a  system  of  integral  equations,  we  first  solve  for  4/””,  (p'/" 
in  the  linear  system  involving  obtaining  a  new  problem  which  does  not  contain 

dependencies  among  the  variables.  We  reduce  this  new  problem  to  a  first  order  system,  and 
then  use  the  results  of  Theorem  2.10  to  express  $  by  the  formula 

$(i)=  4'(a:).r(x), 

with  'J'  :  [a,c]  — *  given  by  the  formula 


(223) 
(22-1 ) 
(223) 


3'(x)  = 


(220) 


with  4'a,  4'b  ■  [-2  5  2]  L(R^^^)  defined  by  the  formulae 


'I’Aix) 


V>b(x) 


(22?) 


(22S) 


and  with  F  :  [a,c]  — *  R"  the  solution  to  an  equation  of  the  form  (40),  subject  to  boundary 
conditions  of  the  form  (41).  The  results  of  Theorem  2.4  are  then  applied  to  reduce  this  new 
problem  with  inhomogeneous  boundary  conditions  to  a  problem  which  boundary  conditions  are 
homogeneous.  For  these  experiments,  the  background  Green’s  function  is  chosen  to  correspond 
to  the  equation 

$'(.t)  =  0,  (223) 


subject  to  boundary  conditions  (20G). 
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